Regresie liniara clasica
Cozma Laura-Elena (Gr 405)
Iamandii Ana-Maria (Gr 405)
Manolache Andrei (Gr 405)
22 January 2022
1. Regresia Liniara Simplă
1.1. Definirea regresiei liniare simple
Regresia liniara simpla este o metoda statistica prin care se descrie relatia dintre doua variabile:
- Variabila independenta, de obicei notata cu X, numita variabila explicativa sau predictor
- Variabila continua dependenta rezultata, notata cu Y, numita variabila raspuns.
1.1.1 Un exemplu de relație liniară
Pentru a exemplifica o relatie liniara, se va folosi data frame-ul Student Survey Data, din cadrul pachetul MASS. Componentele utilizate in continuare in cadrul exemplului sunt:
- Wr.Hnd, intinderea palmei
- Height, inaltimea persoanei in cm
library(MASS)
plot(survey$Height~survey$Wr.Hnd,xlab="Intinderea palmei (cm)", ylab="Inaltime (cm)",
sub = "Graficul relatiei dintre inaltimea si intinderea palmei studentilor analizati in setul de date")
Dupa cum poate fi observat din graficul de mai sus, exista o asociere pozitiva dintre intinderea palmei si inaltime. Aceasta asociere se poate observa si din calculul coeficientului de corelatie.
cor(survey$Wr.Hnd,survey$Height,use="complete.obs")## [1] 0.60099091.2. Concepte generale
Scopul modelului de regresie liniara este de a gasi o functie pentru a estima valoarea medie a variabilei rezultat, fiind data valoarea predictorului.
In cazul setului de date amintit anterior, se poate considera regresia liniara pentru a determina urmatorul rezultat: ce inaltime ne asteptam sa aiba un student care are intinderea mainii de 15 cm?
1.2.1 Definirea modelului
Modelul regresiei liniare simple defineste valoarea raspunsului conditionat de valoarea variabilei explicative cu ajutorul urmatoarei formule:
\[ Y|X = \beta_0 + \beta_1X + \epsilon \]
Parametrii modelului
- X - predictorul sau variabila explicativa
- Y - variabila dependenta
- \(\beta_0\) - interceptorul - valoarea rezultata cand predictorul e 0
- \(\beta_1\) - panta - schimbarea pe care o are valoarea medie rezultata odata cu cresterea predictorului
- \(\epsilon\) - reziduul sau eroarea estimarii
Regresia liniară calculeaza dreapta cu cea mai bună potrivire pentru datele de antrenare, gasind coeficientii de regresie care minimizează eroarea totală \(\epsilon\) a modelului.
Considerand parametrii estimati \(\hat{\beta_0}\) si \(\hat{\beta_1}\), modelul determina valoarea medie a raspunsului, notata cu \(\hat{y}\), pentru valoarea predictor x, dupa formula:
\[ \hat{y} = \hat{\beta_0} + \hat{\beta_1}x \]
1.2.2 Regresia liniara simpla folosind R
R pune la dispozitie comanda lm care realizeaza estimarea parametrilor:
survfit <- lm(Height~Wr.Hnd,data=survey)
survfit##
## Call:
## lm(formula = Height ~ Wr.Hnd, data = survey)
##
## Coefficients:
## (Intercept) Wr.Hnd
## 113.954 3.117Astfel, putem observa ca modelul estimat este:
\[ \hat{y} = 113.954 + 3.117x \] Conform ecuatiei, putem interpreta ca la fiecare crestere cu 1 cm a intinderii palmei, inaltimea persoanei va creste cu 3,117 cm.
Anterior, am definit reziduul ca fiind eroarea estimarii. Aceasta poate fi vazuta pe grafic ca segmentul trasat intre valoarea estimata de regresia liniara si valoarea reala. Modelul de regresie liniara simpla vrea sa minimizeze aceste erori, ea fiind dreapta ce se va afla cel mai aproape de toate observatiile.
In continuare vom reprezenta graficul observatiilor din setul de date, impreuna cu dreapta determinata de ecuatia modelului de regresie. De asemenea, vom trasa si segmentul dintre valorile reale si valoarile estimate a doua observatii, pentru a ilustra eroarea estimarii.
plot(survey$Height~survey$Wr.Hnd,xlab="Intinderea mainii (cm)", ylab="Inaltime (cm)",
sub = "Linia regresiei liniare simple, erorile de estimare si graficul datelor din survey ")
abline(survfit,lwd=2)
obsA <- c(survey$Wr.Hnd[197],survey$Height[197])
obsB <- c(survey$Wr.Hnd[154],survey$Height[154])
mycoefs <- coef(survfit)
resid <- mycoefs[1]
fitted <- mycoefs[2]
segments(x0 = c(obsA[1], obsB[1]), y0 = resid + fitted * c(obsA[1], obsB[1]), x1 = c(obsA[1], obsB[1]), y1 = c(obsA[2], obsB[2]), lty=2)
1.3. Inferenta statistica
In cadrul regresiei liniare simple, trebuie sa ne intrebam daca exista o dovada care sa sustina faptul ca o schimbare in variabila explicativa afecteaza raspunsul mediu. In continuare vom vedea mai multe metode de testare a acestei ipoteze.
1.3.1 Metoda summary din R
Inferenta bazata pe model este realizata direct in R cu ajutorul metodei summary. Vom aplica metoda pe modelul nostru:
summary(survfit)##
## Call:
## lm(formula = Height ~ Wr.Hnd, data = survey)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.7276 -5.0706 -0.8269 4.9473 25.8704
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 113.9536 5.4416 20.94 <2e-16 ***
## Wr.Hnd 3.1166 0.2888 10.79 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.909 on 206 degrees of freedom
## (29 observations deleted due to missingness)
## Multiple R-squared: 0.3612, Adjusted R-squared: 0.3581
## F-statistic: 116.5 on 1 and 206 DF, p-value: < 2.2e-16Reziduul
Reziduul reprezinta diferenta dintre valoarea actuala si cea estimata. Valorea medianei ne spune daca daca exista o simetrie in reziduurile rezultate.
Coeficientii
Prima coloana din sectiunea coefficients contine valorile estimate pentru interceptor si panta.
Urmeaza Std. Error, deviatia standard. Aceasta estimeaza deviatia standard a coeficientilor - cat de multa incertitudine exista in estimarea coeficientilor. Valorile pot fi folosite pentru a afla intervale de incredere.
Valoarea t este rezultata prin calcularea raportului dintre coeficient si deviatia standard. In general, vrem ca valoarea t sa fie mare, deoarece indica faptul ca deviatia standard este mica fata de coeficient. Astfel, cu cat este mai mare valoarea t, cu atat putem fi mai siguri ca valoarea coeficientului nu este 0.
Valoarea p ne ajuta sa intelegem cat de semnificativ este coeficientul nostru in model. De obicei, valorile sub 0.05 sunt considerate semnificative. In cazul in care acesta este semnificativ, inseamna ca ar adauga valoare in estimarea rezultatului. De asemenea, numarul de stelute din dreptul valorii p ne indica cat de semnificativ este un parametru.
Eroarea reziduala standard - residual standard error
Aceasta masura ne spune cat de bine a estimat modelul datele. Din rezultat, putem deduce faptul ca valorile reale ale inaltimii sunt, in medie, cu 7.9 cm mai departe fata de valorile rezultate.
Multiple R-squared
Valoarea lui multiple R-squared descrie ce procent din variatia raspunsului este influentat de predictor. Aceasta este o metoda prin care se determina cat de bine se potriveste modelul cu datele noastre. In cazul regresiei liniare simple, multiple R-squared poate fi calculat si ca valoarea corelatiei la patrat.
In exemplul nostru, 36% din variatia din cadrul inaltimii este determinata de intinderea mainii.
F-statistic si valoarea p
F-statistic si valoarea lui p ne poate ajuta si pentru a confirma sau infirma ipoteza nula. Ipoteza nula afirma ca nu exista o relatie intre variabila dependenta si cea independenta. Vom defini ipotezele:
\[ H_0: \beta_j = 0 \\ H_A: \beta_j \neq 0 \]
Noi dorim sa stim daca exista un efect dat de predictor, adica panta este diferita de 0. Este cunoscut faptul ca daca valoarea p este cu mult mai mica ca 0, atunci exista o dovada puternica asupra faptului ca ipoteza nula este falsa, astfel putem concluziona faptul ca exista o relatie intre variabila independenta si cea dependenta.
In exemplul nostru, valoarea p este \(2 * 10^{-16}\), deci ipoteza nula este infirmata: o crestere a intinderii mainii este asociata cu o crestere in inaltime a populatiei care este studiata. De asemenea, pentru modele mici, o valoare F-statistic mare ne indica faptul ca ipoteza nula poate fi respinsa.
1.4. Predictia valorilor
Cand dorim sa calculam rezultatul unei predictii, inlocuim in ecuatie valoarea lui x de care suntem interesati, insa rezultatul final va fi supus unei variatii,
1.4.1 Intervale de incredere si intervale de predictie
Intervalele de confidenta sunt folosite pentru a descrie variabilitatea mediei raspunsului, in timp ce intervalele de predictie sunt folosite pentru a oferi intervalul posibil de valori pe care o observatie individuala o poate avea, fiind dat x.
In continuare vom calcula intervalele de incredere si de predictie pentru un handspan de 14.5 cm si de 24 cm cu un nivel de incredere de 95%, folosind functia predict din R:
xvals <- data.frame(Wr.Hnd=c(14.5,24))
mypred.ci <- predict(survfit,newdata=xvals,interval="confidence",level=0.95)
mypred.ci## fit lwr upr
## 1 159.1446 156.4956 161.7936
## 2 188.7524 185.5726 191.9323mypred.pi <- predict(survfit,newdata=xvals,interval="prediction",level=0.95)
mypred.pi## fit lwr upr
## 1 159.1446 143.3286 174.9605
## 2 188.7524 172.8390 204.6659plot(survey$Height~survey$Wr.Hnd,xlim=c(13,24),ylim=c(140,205),
xlab="Writing handspan (cm)",ylab="Height (cm)")
abline(survfit,lwd=2)
points(xvals[,1],mypred.ci[,1],pch=8)
segments(x0=c(14.5,24),y0=c(mypred.pi[1,2], mypred.pi[2,2]), x1=c(14.5,24), y1=c(mypred.pi[1,3],mypred.pi[2,3]),col="gray",lwd=3)
segments(x0=c(14.5,24),y0=c(mypred.ci[1,2], mypred.ci[2,2]), x1=c(14.5,24), y1=c(mypred.ci[1,3], mypred.ci[2,3]),lwd=2)
xseq <- data.frame(Wr.Hnd=seq(12,25,length=100))
ci.band <- predict(survfit,newdata=xseq,interval="confidence",level=0.95)
pi.band <- predict(survfit,newdata=xseq,interval="prediction",level=0.95)
lines(xseq[,1],ci.band[,2],lty=2)
lines(xseq[,1],ci.band[,3],lty=2)
lines(xseq[,1],pi.band[,2],lty=2,col="gray")
lines(xseq[,1],pi.band[,3],lty=2,col="gray")
legend("topleft",legend=c("Fit","95% CI","95% PI"),lty=c(1,2,2),
col=c("black","black","gray"),lwd=c(2,1,1))
1.4.2 Interpolarea si Extrapolarea
Interpolarea si extrapolarea sunt doi termeni ce descriu natura unei predictii.
O predictie este referita ca fiind interpolare cand valoarea x a predictorului se afla in intervalul datelor observate la antrenare. In momentul in care x este in exteriorul intervalului, are loc o extrapolare.
In exemplul analizat anterior, cazul x = 14.5 cm este o interpolare, iar x = 24 cm este o extrapolare.
In general, este preferabil sa fie utilizata interpolarea, deoarece se alege o valoare din intervalul care a fost analizat, sau extrapolarea in cazul in care x este in vecinatatea intervalului.
Un exemplu de extrapolare care ofera rezultate nerealiste este x = 0. Astfel, un om cu intinderea mainii de 0 cm ar avea o inaltime egala cu 114 cm, rezultatele fiind improbabile. Un alt exemplu aberant ar fi pentru o intindere a mainii de 55 cm. Vom calcula rezultatul folosind functia predict.
predict(survfit,newdata=data.frame(Wr.Hnd=55),interval="confidence",
level=0.95)## fit lwr upr
## 1 285.3675 264.6992 306.0359Se observa ca valoarea prezisa este de 285 cm, ceea ce nu are sens in contextul problemei analizate.
Exercitiul 20.1
a.
pred_a.ci <- predict(survfit, newdata = data.frame(Wr.Hnd = c(12, 15.3, 17, 19.9)), interval = "confidence", level = 0.99)
pred_a.ci## fit lwr upr
## 1 151.3530 146.0901 156.6159
## 2 161.6379 158.6827 164.5930
## 3 166.9361 164.9985 168.8737
## 4 175.9743 174.3066 177.6420b.
# ia observatiile incomplete
incomplete.obs <- which(is.na(survey$Height) | is.na(survey$Wr.Hnd))
rho.xy <- cor(survey$Wr.Hnd, survey$Height, use = "complete.obs")
# determina coeficientii cu ajutorul formulelor
b1 <- rho.xy * (sd(survey$Height[-incomplete.obs]) / sd(survey$Wr.Hnd[-incomplete.obs]))
b0 <- mean(survey$Height[-incomplete.obs])-b1*mean(survey$Wr.Hnd[-incomplete.obs])
b0## [1] 113.9536b1## [1] 3.116617# afiseaza coeficientii rezultati folosind survfit
survfit##
## Call:
## lm(formula = Height ~ Wr.Hnd, data = survey)
##
## Coefficients:
## (Intercept) Wr.Hnd
## 113.954 3.117c. i.
Ecuatia este: y = 177.857 - 0.072x
survfit.pulse <- lm(Height~Pulse, data=survey)
survfit.pulse##
## Call:
## lm(formula = Height ~ Pulse, data = survey)
##
## Coefficients:
## (Intercept) Pulse
## 177.85708 -0.07225plot(survey$Height~survey$Pulse, xlab="Pulsul (bpm)",ylab="Inaltimea (cm)")
abline(survfit.pulse,lwd=2)
summary(survfit.pulse)##
## Call:
## lm(formula = Height ~ Pulse, data = survey)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.3543 -7.2019 -0.9439 7.2622 26.1168
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 177.85708 4.93485 36.041 <2e-16 ***
## Pulse -0.07225 0.06598 -1.095 0.275
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.884 on 169 degrees of freedom
## (66 observations deleted due to missingness)
## Multiple R-squared: 0.007046, Adjusted R-squared: 0.001171
## F-statistic: 1.199 on 1 and 169 DF, p-value: 0.275confint(survfit.pulse, level = 0.9)## 5 % 95 %
## (Intercept) 169.6952304 186.01892281
## Pulse -0.1813757 0.03687061Conform rezultatelor, pentru fiecare bataie pe minut in plus, media inaltimii studentului scade cu 0.072.
Din calcularea intervalelor de incredere, se observa ca intervalul de incredere contine valoarea 0, deci panta ar putea sa fie 0. De asemenea, valoarea p pentru panta este 0.275, nu suficient de mica incat sa fie considerata semnificativa. Nu exista dovezi care ar putea respinge H0, astfel ca nu poate fi concluzionat faptul ca pulsul afecteaza inaltimea unui student. Aceasta ipoteza este sustinuta si printr-un F-statistic mic de 1.199, dar si prin multiple R-squared cu valoarea mica de 0.007, astfel ca numai 0.7% din inaltime e influentata de puls.
# definim un interval de interes
x_seq <- data.frame(Pulse=seq(45,120,length=100))
plot(survey$Height~survey$Pulse, xlab="Pulsul (bpm)",ylab="Inaltimea (cm)")
pred_iii.ci <- predict(survfit.pulse, newdata = x_seq, interval = "confidence", level = 0.9)
pred_iii.pi <- predict(survfit.pulse, newdata = x_seq, interval = "prediction", level = 0.9)
lines(x_seq[,1], pred_iii.ci[,2], lty=2)
lines(x_seq[,1], pred_iii.ci[,3], lty=2)
lines(x_seq[,1], pred_iii.pi[,2], lty=2, col="gray")
lines(x_seq[,1], pred_iii.pi[,3], lty=2, col="gray")
legend("topright",legend=c("CI 90%","PI 90%"),lty=2,col=c("black","gray"))
# creeaza vectorul incomplete
incomplete.obs <- which(is.na(survey$Height) | is.na(survey$Pulse))
plot(survey$Height~survey$Pulse, xlab="Pulsul (bpm)",ylab="Inaltimea (cm)")
m <- mean(survey$Height[-incomplete.obs])
abline(h = m, col='brown', lty=3, lwd=3)
pred_iii.ci <- predict(survfit.pulse, newdata = x_seq, interval = "confidence", level = 0.9)
pred_iii.pi <- predict(survfit.pulse, newdata = x_seq, interval = "prediction", level = 0.9)
lines(x_seq[,1], pred_iii.ci[,2], lty=2)
lines(x_seq[,1], pred_iii.ci[,3], lty=2)
lines(x_seq[,1], pred_iii.pi[,2], lty=2, col="gray")
lines(x_seq[,1], pred_iii.pi[,3], lty=2, col="gray")
Dreapta in interiorul intervalului de incredere, ceea ce sustine faptul ca panta ar avea valoarea 0, deci pulsul nu influenteaza inaltimea unui student.
d.
?mtcars## starting httpd help server ... doneplot(mtcars$mpg~mtcars$wt, xlab="Greutate", ylab="MPG")
e.
cars.fit<- lm(mpg~wt,data=mtcars)
cars.fit##
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
##
## Coefficients:
## (Intercept) wt
## 37.285 -5.344plot(mtcars$mpg~mtcars$wt, xlab="Greutate", ylab="MPG")
abline(cars.fit, lwd = 2)
f.
summary(cars.fit)##
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5432 -2.3647 -0.1252 1.4096 6.8727
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.2851 1.8776 19.858 < 2e-16 ***
## wt -5.3445 0.5591 -9.559 1.29e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.046 on 30 degrees of freedom
## Multiple R-squared: 0.7528, Adjusted R-squared: 0.7446
## F-statistic: 91.38 on 1 and 30 DF, p-value: 1.294e-10Ecuatia este: y = 37.285 - 5.344x, unde y reprezinta MPG mediu si x este greutatea.
La fiecare crestere a greutatii cu o unitate (1000 lbs), mpg scade cu -5.344. Valoarea lui p este foarte mica (1.29e-10), deci ipoteza nula este infirmata, ceea ce arata ca exista o dovada statistica puternica ca MPG se schimba odata cu greutatea.
g.
predict(cars.fit, newdata = data.frame(wt = 6), interval = "prediction", level = 0.95)## fit lwr upr
## 1 5.218297 -1.85279 12.28938Modelul nu prezice cu exactitate observatiile, deoarece a fost folosita o extrapolare cu mult in exteriorul intervalului de date observate. Astfel, PI asociat are o limita inferioara negativa, care nu are sens in contextul MPG.
1.5. Predictori categorici
In sectiunea anterioara a fost analizat modelul regresiei liniare simple asupra unei variabile explicative continue. In continuare se va analiza cazul in care variabila explicativa este categorica - o variabila care poate lua un numar fix de posibile valori, fiecare reprezentand un grup sau o categorie.
1.5.1 Variabile binare: k = 2
Vom trata cazul variabilelor binare, in care predictorul va avea doua categorii cu valorile 0, respectiv 1.
Ecuatia regresiei liniare simple in cazul in care variabila predictor e categorica se pastreaza, aceasta fiind:
\[Y|X=\beta_0 + \beta_1X + \epsilon\]
Insa semnificatiile coeficientilor ecuatiilor se vor schimba, devenind:
- \(\beta_0\) va avea semnificatia de baseline sau referinta pentru cazul in care \(X = 0\), intrucat Y devine \(Y = \beta_0 + \epsilon\)
- \(\beta_1\) va fi efectul aditiv in cazul in care \(X = 1\), deoarece \(Y = \beta_0 + \beta_1 + \epsilon\)
Analiza grafica
Pentru a analiza modelul de regresie liniara simpla asupra variabilelor categorice, se vor folosi variabilele Sex - de tip factor, cu valorile Female si Male - si Height din data frame-ul Survey. Vrem sa determinam daca exista o dovada statistica conforma careia inaltimea unui student este influentata de sexul acestuia.
Vom construi un plot ce contine doua boxplot-uri cu cele doua categorii, impreuna cu punctele determinate de datele de intrare si cu inaltimile medii.
plot(survey$Height~survey$Sex, xlab = "Sex", ylab="Inaltime (cm)", sub="Graficul relatiei dintre inaltime studentilor si categoria in care se afla: Female sau Male")
points(survey$Height~as.numeric(survey$Sex),cex=0.5)
means.sex <- tapply(survey$Height,INDEX=survey$Sex,FUN=mean,na.rm=TRUE)
points(1:2,means.sex,pch=4,cex=2)
In urma plotului, se observa ca barbatii tind sa fie mai inalti decat femeile, dar trebuie gasite dovezi statistice care sa sustina aceasta afirmatie.
Regresie liniara cu predictor categoric folosind R
Vom folosi functia lm pentru a estima coeficientii:
survfit2 <- lm(Height~Sex,data=survey)
summary(survfit2)##
## Call:
## lm(formula = Height ~ Sex, data = survey)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23.886 -5.667 1.174 4.358 21.174
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 165.687 0.730 226.98 <2e-16 ***
## SexMale 13.139 1.022 12.85 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.372 on 206 degrees of freedom
## (29 observations deleted due to missingness)
## Multiple R-squared: 0.4449, Adjusted R-squared: 0.4422
## F-statistic: 165.1 on 1 and 206 DF, p-value: < 2.2e-16Spre deosebire de estimarea regresiei cu variabila independenta continua, \(\beta_0\), interceptorul, va reprezenta acum valoarea inaltimii medii estimate in cazul in care studentul are genul feminin (x = 0).
\(\beta_1\) este raportat ca SexMale, astfel ca 13.139 este diferenta in inaltime care este adaugata in cazul in care variabila x defineste genul masculin (x = 1).
Valoarea lui \(\beta_1\) va determina dovada statistica prin care media variabilei de raspuns este afectata de variabila explicativa, astfel ca cu cat \(\beta_1\) e mai mare, cu atat rezultatele vor mai diferite pentru cele doua cazuri.
Predictia modelului
In urma rezultatului dat de functia lm, a rezultat urmatoarea ecuatie ce defineste modelul:
\[ \hat{y} = 165.687 + 13.139x \]
Vom defini un vector cu valorile Male si Female pentru a realiza predictia rezultatelor.
extra.obs <- factor(c("Female","Male","Male","Male","Female"))
extra.obs## [1] Female Male Male Male Female
## Levels: Female Malepredict(survfit2,newdata=data.frame(Sex=extra.obs),interval="confidence",
level=0.9)## fit lwr upr
## 1 165.6867 164.4806 166.8928
## 2 178.8260 177.6429 180.0092
## 3 178.8260 177.6429 180.0092
## 4 178.8260 177.6429 180.0092
## 5 165.6867 164.4806 166.8928Se observa ca valorile predictiilor difera numai intre cele doua categorii, fiind 165.68 pentru Female si 178.82 pentru Male, intrucat ecuatiile rezolvate sunt:
- \(\hat{y} = 165.68\), pentru Female, x = 0
- \(\hat{y} = 165.68 + 13.139\), pentru Male, x = 1
1.5.2 Variabile cu nivel multiplu: k > 2
Vom trata cazul in care variabila predictor are mai mult de 2 niveluri. Se va utiliza reprezentarea:
\[ X_{(1)} = 0,1; X_{(2)} = 0,1; X_{(3)} = 0,1; ...; X_{(k)} = 0,1 \] Aceasta procedura se numeste dummy coding.
Astfel, daca un individ a ales \(X = 2\) pentru variabila categorica, atunci \(X_{(2)} = 1\) si toate celelalte valori sunt 0.
Ideea acesteia se bazeaza pe faptul ca intre variabilele categorice nu pot fi, in mod general, create legaturi numerice ca in cazul celor continue. Astfel, in cadrul modelului a fost abordat cazul absentei sau prezentei binare.
Formula modelului de regresie liniara pentru un predictor categoric
Formula modelului de regresie liniara simpla pentru o variabila independenta categorica avand k > 2 este:
\[ \hat{y} = \hat{\beta_0} + \hat{\beta_1}X_{(2)} + \hat{\beta_2}X_{(3)} + . . . + \hat{\beta}_{k-1}X_{(k)} \] Prima categorie a predictorului este considerata nivel de referinta. Drept exemplu, in cazul in care categoria aleasa este 3, \(X_{(3)} = 1\), toate celelalte sunt 0, iar ecuatia devine: \(\hat{y} = \hat{\beta_0} + \hat{\beta_2}\)
Regresia liniara pentru variabile cu nivel multiplu folosind R
Se va analiza daca inaltimea unui individ din setul de date este influentata de frecventa fumatului, data prin intermediul variabilei Smoke.
boxplot(Height~Smoke,data=survey)
points(1:4,tapply(survey$Height,survey$Smoke,mean,na.rm=TRUE),pch=4)
Calculam predictia coeficientilor cu ajutorul functiei lm:
survfit3 <- lm(Height~Smoke,data=survey)
summary(survfit3)##
## Call:
## lm(formula = Height ~ Smoke, data = survey)
##
## Residuals:
## Min 1Q Median 3Q Max
## -25.02 -6.82 -1.64 8.18 28.18
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 173.7720 3.1028 56.005 <2e-16 ***
## SmokeNever -1.9520 3.1933 -0.611 0.542
## SmokeOccas -0.7433 3.9553 -0.188 0.851
## SmokeRegul 3.6451 4.0625 0.897 0.371
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.812 on 205 degrees of freedom
## (28 observations deleted due to missingness)
## Multiple R-squared: 0.02153, Adjusted R-squared: 0.007214
## F-statistic: 1.504 on 3 and 205 DF, p-value: 0.2147Se observa ca Heavy a fost considerata categoria de referinta, coeficientul acesteia fiind reprezentata de interceptor.
Predictia valorilor
Se va realiza cu ajutorul functiei predict:
one.of.each <- factor(levels(survey$Smoke))
one.of.each## [1] Heavy Never Occas Regul
## Levels: Heavy Never Occas Regulpredict(survfit3,newdata=data.frame(Smoke=one.of.each),
interval="confidence",level=0.95)## fit lwr upr
## 1 173.7720 167.6545 179.8895
## 2 171.8200 170.3319 173.3081
## 3 173.0287 168.1924 177.8651
## 4 177.4171 172.2469 182.5874Rezultatul dat de metoda summary arata ca nicio valoare p corespunzatoare coeficientilor nu este suficient de semnificativa. Se poate concluziona ca inaltimea medie a unui student nu este influentata de frecventa fumatului.
Aceasta ipoteza este intarita si de valoarea mica a multiple r-squared, ce indica o variatie mica a inaltimii influentata de frecventa fumatului.
Desi coeficientul ${_0} (173.772) este destul de mare, nu poate sugera decat faptul ca inaltimea cel mai probabil nu poate fi 0.
1.5.3 Tratarea variabilelor categorice ca variabile numerice
In cazul in care datele categorice ce trebuie analizate nu sunt sub forma unui factor in cadrul setului de date si sunt sub forma numerica, lm va trata variabila ca un predictor continuu, schimband per unitate raspunsul mediu. Aceasta abordare nu este potrivita daca variabila explanatorie are rolul de a reprezenta grupuri distincte.
Desi aceasta metoda poate părea inadecvata daca variabila explicativa initiala se presupune a fi alcatuita din grupuri distince, in unele situatii, cand variabila poate fi tratata în mod natural ca fiind discreta, acest tratament este valabil statistic si ajută la interpretare.
Analiza grafica
Se analizeaza setul de date mtcars, unde variabila de raspuns e kilometrajul mpg, iar predictorul este numarul de cilindri cly. Desi cly este o variabila categorica, ea ia valori numerice si nefiind un vector factor, lm o va trata ca fiind o variabila continua continua. Conform parametrilor determinati, se considera ca la fiecare cilindru adaugat, kilometrajul scade cu 2.88 MPG in medie.
carfit <- lm(mpg~cyl,data=mtcars)
summary(carfit)##
## Call:
## lm(formula = mpg ~ cyl, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.9814 -2.1185 0.2217 1.0717 7.5186
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.8846 2.0738 18.27 < 2e-16 ***
## cyl -2.8758 0.3224 -8.92 6.11e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.206 on 30 degrees of freedom
## Multiple R-squared: 0.7262, Adjusted R-squared: 0.7171
## F-statistic: 79.56 on 1 and 30 DF, p-value: 6.113e-10Conform rezultatului din *summary, ecuatia modelului este: \[ \hat{y} = 37.88 - 2.88x \]
Astfel, la fiecare cilindru adaugat, kilometrajul scade cu 2.88 MPG in medie.
boxplot(mtcars$mpg~mtcars$cyl,xlab="Cilindri",ylab="MPG")
plot(mtcars$mpg~mtcars$cyl,xlab="Cilindri",ylab="MPG")
abline(carfit,lwd=2)
Utilizarea predictorilor categorici ca variabile continue prezinta o serie de avantaje:
- Poate fi folosita la interpolare - se poate evalua valoarea medie MPG pentru o masina cu 5 cilindri.
- Necesita aflarea unui numar mai mic de parametri - doar panta in locul a k-1 coeficienti.
Totusi, exista si dezavantaje:
- Nu se primesc informatiile la nivel de grup
- Raspunsul mediu corespunzator fiecarei categorii poate sa nu fie reprezentat bine cu abordarea continua
Echivalenta cu ANOVA
Ceea ce realizeaza regresia liniara simpla cu un singur predictor categoric este echivalenta cu metoda ANOVA unidirectionala: compararea a mai mult de doua valori medii si determinarea unei unei dovezi statistice conform careia cel putin o medie se distinge fata de celelalte.
Procedura ANOVA determina valoarea p globala care care cuantifica nivelul dovezilor statistice asupra ipotezei nule. Aceasta valoare rezulta si la finalul rezultatului metodei summary aplicata asupra obiectului lm.
Exercitiul 20.2
a.
table(survey$Exer)##
## Freq None Some
## 115 24 98boxplot(survey$Height~survey$Exer, xlab="Exercitii fizice",
ylab="Inaltime(cm)")
b.
exer.fit <- lm(Height~Exer, data = survey)
summary(exer.fit)##
## Call:
## lm(formula = Height ~ Exer, data = survey)
##
## Residuals:
## Min 1Q Median 3Q Max
## -24.607 -6.397 -1.607 6.103 25.393
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 174.6067 0.9396 185.836 < 2e-16 ***
## ExerNone -5.5787 2.3489 -2.375 0.01847 *
## ExerSome -4.2098 1.4094 -2.987 0.00316 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.628 on 206 degrees of freedom
## (28 observations deleted due to missingness)
## Multiple R-squared: 0.05333, Adjusted R-squared: 0.04414
## F-statistic: 5.802 on 2 and 206 DF, p-value: 0.003536Nivelul de referinta considerat implicit este Freq, primul in ordine alfabetica.
c.
Conform predictiei, coeficientii sunt: Freq: 174.606 None: -5.5787
Some: -4.2098 Deoarece Freq este considerat nivelul de referinta, modelul va prezice efectul pe care il va produce None sau Some asupra inaltimii. Avand in vedere ca None si Some au coeficienti negativi, din model se deduce ca persoanele care exerseaza mai putin decat frecvent, vor fi in medie mai scunde, intrucat din valoarea de referinta se va scadea. Deoarece -5.57 este cel mai mic numar negativ dintre coeficienti, inseamna ca persoanele care nu exerseaza deloc vor avea inaltimea cea mai mica (174.606 - 5.5787 = 169.02 cm)
Analizand valorile p corespunzatoare coeficientilor, atat pentru ExerNone cat si pentru ExerSome sunt suficient de mici incat coeficientii sa fie considerati statistic semnificativ diferiti de 0.
De asemenea, valoarea p globala este 0.0035, suficient de mica cat poata fi afirmat faptul ca exercitiile fizice influenteaza inaltimea.
d.
lev <- factor(levels(survey$Exer))
lev## [1] Freq None Some
## Levels: Freq None Somepredict(exer.fit,newdata = data.frame(Exer = lev), interval = "prediction", level = 0.95)## fit lwr upr
## 1 174.6067 155.5349 193.6784
## 2 169.0280 149.5777 188.4783
## 3 170.3969 151.3027 189.4911e.
anova <- aov(Height~Exer, data = survey)
summary(anova)## Df Sum Sq Mean Sq F value Pr(>F)
## Exer 2 1076 537.8 5.802 0.00354 **
## Residuals 206 19095 92.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 28 observations deleted due to missingnessIn urma analizei ANOVA se obtine aceeasi valoare p globala egala cu 0.0035, astfel sustinandu-se ideea ca exercitiile fizice afecteaza inaltimea unui individ din grupul studiat.
f.
exer.new_order <- relevel(survey$Exer,ref="None")
levels(exer.new_order)## [1] "None" "Freq" "Some"anova <- aov(Height~Exer, data = survey)
summary(anova)## Df Sum Sq Mean Sq F value Pr(>F)
## Exer 2 1076 537.8 5.802 0.00354 **
## Residuals 206 19095 92.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 28 observations deleted due to missingnessValoarea p globala determinata de ANOVA este aceeasi. Acesta este si rezultatul asteptat, intrucat o schimbare in ordinea factorilor nu trebuie sa afecteze rezultatele prezise.
g.
car_qsec.fit <- lm(qsec~gear, data = mtcars)
summary(car_qsec.fit)##
## Call:
## lm(formula = qsec ~ gear, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.7929 -1.1604 -0.3278 1.2122 5.2122
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.7482 1.6239 12.161 4e-13 ***
## gear -0.5151 0.4321 -1.192 0.243
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.775 on 30 degrees of freedom
## Multiple R-squared: 0.04523, Adjusted R-squared: 0.01341
## F-statistic: 1.421 on 1 and 30 DF, p-value: 0.2425Desi gear este o variabila categorica, avand valori numerice care nu sunt stocate sub forma unui factor, aceasta va fi tratata ca o variabila continua. Astfel, pentru fiecare crestere a vitezei cu o unitate, exista o scadere a timpului cu 0.51 secunde.
Insa valoarea p globala nu este suficient de mica (0.242), astfel ca nu exista dovezi statistice care sa sustina faptul ca schimbarea vitezei determina o schimbare a timpului.
h.
car_qsec.fit_cat <- lm(qsec~factor(gear),data=mtcars)
summary(car_qsec.fit_cat)##
## Call:
## lm(formula = qsec ~ factor(gear), data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5050 -0.6667 -0.2060 0.6125 3.9350
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.6920 0.3691 47.928 < 2e-16 ***
## factor(gear)4 1.2730 0.5537 2.299 0.02890 *
## factor(gear)5 -2.0520 0.7383 -2.779 0.00946 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.43 on 29 degrees of freedom
## Multiple R-squared: 0.4012, Adjusted R-squared: 0.3599
## F-statistic: 9.715 on 2 and 29 DF, p-value: 0.0005897In cazul in care gear e folosit ca variabila categorica, valoarea p globala este mult mai mica, 0.00058, aducand dovada statistica ca gear influenteaza qsec. De asemenea, si valorile p corespunzatoare nivelurilor 3 (cel dat de interceptor), 4 si 5 sunt suficient de mici cat sa poata fi considerati coeficientii ca aducand valoare rezultatului final.
Astfel, pentru gear = 3, qsec = 17.69, gear = 4 aduce o valoare aditiva de 1.27 secunde, iar gear = 5 scade qsec cu 2.05 secunde.
i.
boxplot(mtcars$qsec~mtcars$gear,xlab="Gears",ylab="Qsec")
plot(mtcars$qsec~mtcars$gear,xlab="Gears",ylab="Qsec")
abline(car_qsec.fit,lwd=2)
Cele doua ploturi ilustreaza faptul ca regresia liniara trebuie sa trateze gears ca pe o variabila categorica si nu una continua, intrucat relatia nu poate fi explicata cu ajutorul unei drepte.
Astfel, abordarea folosind variabile continue nu va reprezenta in mod corect schimbarea pe care o are gears asupra qsec.
2. Regresia Liniara Multipla
Regresia liniara multipla este o generalizare directa a modelelor cu un singur predictor discutate pana acum. Aceasta permite modelarea variabilei de raspuns in functie de mai multi predictori, astfel incat sa se poate observa efectul variabilelor explicative asupra variabilei de raspuns.
2.1. Terminologie
Variabila lurking: influenteaza raspunsul, alt predictor sau ambele, dar nu este inclusa intr-un model de predictie.
Confounding: Prezenta unei variabile lurking poate duce la concluzii false despre relatiile de cauzalitate dintre raspuns si alt predictor sau poate masca o asociere cauza-efect adevarata. Aceasta eroare se numeste confounding.
2.2. Formula
Vrem sa determinam valoarea unei variabile de raspuns continue Y fiind date valorile a p>1 variabile explicative independente \[ X_0, X_1, X_2, ... X_p\ \]
\[ Y = \beta_0 + \beta_1\ X_1 + \beta_2\ X_2 + ... + \beta_p\ X_p + \epsilon \]
\[ \epsilon \sim N(0,\sigma^2) \]
unde \[ \beta_0, \beta_1, \beta_2, ... \beta_p\ \] sunt coeficientii regresiei
2.3. Reprezentare cu matrici
\[
Y = X \beta+ \epsilon
\] 
2.4. Implementare si interpretare
Pentru urmatorul exemplu folosim setul de date survey din MASS package.
#library(MASS)
data(survey)
head(survey)## Sex Wr.Hnd NW.Hnd W.Hnd Fold Pulse Clap Exer Smoke Height M.I
## 1 Female 18.5 18.0 Right R on L 92 Left Some Never 173.00 Metric
## 2 Male 19.5 20.5 Left R on L 104 Left None Regul 177.80 Imperial
## 3 Male 18.0 13.3 Right L on R 87 Neither None Occas NA <NA>
## 4 Male 18.8 18.9 Right R on L NA Neither None Never 160.00 Metric
## 5 Male 20.0 20.0 Right Neither 35 Right Some Never 165.00 Metric
## 6 Female 18.0 17.7 Right L on R 64 Right Some Never 172.72 Imperial
## Age
## 1 18.250
## 2 17.583
## 3 16.917
## 4 20.333
## 5 23.667
## 6 21.000Anterior, am observat cateva modele de regresie liniara bazate pe variabila de raspuns inaltimea unui student precum si pe predictorii handspan (intinderea palmei, continuu) sau sex (categoric). Am descoperit ca intinderea palmei era semnificativa din punct de vedere statistic, caci coeficientul estimat arata o crestere medie de aproximativ 3.12 cm pentru fiecare crestere de 1 cm in intinderea palmei. Cand am folosit sexul ca variabila explicativa, modelul a sugerat ca barbatii au in plus 13.14 cm la inaltimea medie comparat cu media femeilor.
Totusi, aceste modele nu ne pot arata efectul comun al sex-ului si intinderii palmei asupra prezicerii inaltimii. Daca le folosim impreuna intr-un model de regresie liniara multipla, putem reduce asa-numita confounding.
survmult <- lm(Height~Wr.Hnd+Sex,data=survey)
summary(survmult)##
## Call:
## lm(formula = Height ~ Wr.Hnd + Sex, data = survey)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.7479 -4.1830 0.7749 4.6665 21.9253
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 137.6870 5.7131 24.100 < 2e-16 ***
## Wr.Hnd 1.5944 0.3229 4.937 1.64e-06 ***
## SexMale 9.4898 1.2287 7.724 5.00e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.987 on 204 degrees of freedom
## (30 observations deleted due to missingness)
## Multiple R-squared: 0.5062, Adjusted R-squared: 0.5014
## F-statistic: 104.6 on 2 and 204 DF, p-value: < 2.2e-16Coeficientul pentru intinderea palmei este acum aprox. 1.59, aproape jumatate din valoarea obtinuta la regresia liniara simpla pentru inaltime. Coeficientul pentru sex s-a redus in comparatie cu regresia liniara simpla si este in continuare semnificativ in prezenta intinderii palmei.
2.5. Interpretarea efectelor marginale
In regresia multipla, estimarea fiecarui predictor ia in considerare efectul tuturor celorlaltor predictori prezenti in model. In cazul nostru, pentru studentii de acelasi sex, cresterea de 1 cm in intinderea palmei conduce la o crestere estimata de 1.5944 cm in inaltimea medie.
Pentru studentii cu intinderea palmei similara, barbatii in medie vor fi cu 9.4898 cm mai inalti decat femeile.
Se observa astfel utilitatea regresiei multiple. In acest exemplu, daca folosim doar regresia liniara simpla, rezultatul este inselator, deoarece o parte din variatia inaltimii este determinata de sex, in timp ce alta este atribuita intinderii palmei.
Modelul poate fi exprimat in urmatorul mod:
Mean height = 137.687 + 1.594 × handspan + 9.49 × sex, unde handspan este exprimat in cm, iar variabila sex are valoarea 1 (barbat), 0 (femeie)
2.6. Vizualizarea modelului
survcoefs <- coef(survmult)
survcoefs## (Intercept) Wr.Hnd SexMale
## 137.686951 1.594446 9.489814as.numeric(survcoefs[1]+survcoefs[3])## [1] 147.1768Poate fi scris ca 2 ecuatii, una pentru femei si una pentru barbati:
- Ecuatia pentru femei:
Mean height = 137.687 + 1.594 × handspan
- Ecuatia pentru barbati:
Mean height = (137.687 + 9.4898) + 1.594 × handspan = 147.177 + 1.594 × handspan
2.7. Plot
survcoefs <- coef(survmult)
plot(survey$Height~survey$Wr.Hnd,
col=c("gray","black")[as.numeric(survey$Sex)],
pch=16,xlab="Writing handspan",ylab="Height")
abline(a=survcoefs[1],b=survcoefs[2],col="gray",lwd=2)
abline(a=survcoefs[1]+survcoefs[3],b=survcoefs[2],col="black",lwd=2)
legend("topleft",legend=levels(survey$Sex),col=c("gray","black"),pch=16)
Initial, realizam un scatterplot cu valorile inaltimii si handspan-ului, impartite pe sex. Apoi, abline adauga linia corespunzatoare femeilor si adauga o a doua linie pentru barbati, bazate pe cele doua ecuatii.
2.8. Predictie pornind de la un model de regresie liniara multipla
Ca exemplu, folosind modelul pentru inaltimea unui student ca functie liniara de handspan si sex, putem estima inaltimea medie a unui student de sex masculin cu o intindere a palmei de 16.5 cm, impreuna cu un interval de incredere.
predict(survmult,newdata=data.frame(Wr.Hnd=16.5,Sex="Male"),
interval="confidence",level=0.95)## fit lwr upr
## 1 173.4851 170.9419 176.0283Rezultatul indica valoarea de 173.48 cm ca inaltime medie a studentului si faptul ca putem fi siguri in proportie de 95% ca valoarea adevarata se afla undeva intre 170.94 si 176.03.
In acelasi mod, inaltimea medie a unei studente cu o intindere a palmei de 13 cm este estimata la 158.42 cm, cu un interval de predictie de 99% intre 139.76 si 177.07 cm.
predict(survmult,newdata=data.frame(Wr.Hnd=13,Sex="Female"),
interval="prediction",level=0.99)## fit lwr upr
## 1 158.4147 139.7611 177.06842.9. Transformarea variabilelor numerice
Uneori, functia liniara definita de ecuatia standard a regresiei poate fi inadecvata in ceea ce priveste relatiile dintre o variabila de raspuns si covariatele selectate. Se poate observa, de exemplu, o curbura in scatterplot intre doua variabile numerice pentru care o linie perfect dreapta nu se potriveste cel mai bine.
Transformarea numerica se refera la aplicarea unei functii matematice asupra valorilor numerice pentru a le redimensiona. Exemplu: radacina patrata a unui numar sau transformarea din Fahrenheit in Celsius.
transformare polinomiala
transformare logaritmica
a) Transformare polinomiala
Pentru a clarifica conceptul curburii polinomiale, consideram intervalul [-4,4] si vectorii urmatori:
x <- seq(-4,4,length=50)
y<-x
y2<-x+x^2
y3<-x+x^2+x^3Acestia au reprezentarile urmatoare:
plot(x,y,type="l", sub="Functie liniara")
plot(x,y2,type="l", sub="Functie patratica")
plot(x,y3,type="l", sub="Functie cubica")
Folosim setul de date mtcars, pentru care consideram variabila disp, care descrie volumul deplasarii motorului in inch cubi si variabila de raspuns mile per galon.
Realizam scatterplot pentru date:
plot(mtcars$disp,mtcars$mpg,xlab="Displacement (cu. in.)",ylab="MPG")
a.1 Regresie liniara simpla
Vrem sa aflam care e cea mai buna reprezentare a relatiei si incepem cu regresia liniara simpla.
car.order1 <- lm(mpg~disp,data=mtcars)
summary(car.order1)##
## Call:
## lm(formula = mpg ~ disp, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.8922 -2.2022 -0.9631 1.6272 7.2305
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 29.599855 1.229720 24.070 < 2e-16 ***
## disp -0.041215 0.004712 -8.747 9.38e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.251 on 30 degrees of freedom
## Multiple R-squared: 0.7183, Adjusted R-squared: 0.709
## F-statistic: 76.51 on 1 and 30 DF, p-value: 9.38e-10Se observa un impact liniar negativ a deplasarii adupra kilometrajului - pentru fiecare inch cubic a volumului deplasarii, media scade cu 0.041 mile per galon.
a.2 Regresie liniara multipla
Adaugam apoi patratul lui disp la model. Observam ca aceasta schimbare e semnificativa din punct de vedere statistic: output-ul care ii corespunde patratului are valoarea p de 0.0031. De asemenea coeficientul de determinare este mai mare decat la prima incercare (0.7927 fata de 0.7183), ceea ce ne indreptateste sa consideram ca al doilea model e mai bun.
car.order2 <- lm(mpg~disp+I(disp^2),data=mtcars)
summary(car.order2)##
## Call:
## lm(formula = mpg ~ disp + I(disp^2), data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.9112 -1.5269 -0.3124 1.3489 5.3946
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.583e+01 2.209e+00 16.221 4.39e-16 ***
## disp -1.053e-01 2.028e-02 -5.192 1.49e-05 ***
## I(disp^2) 1.255e-04 3.891e-05 3.226 0.0031 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.837 on 29 degrees of freedom
## Multiple R-squared: 0.7927, Adjusted R-squared: 0.7784
## F-statistic: 55.46 on 2 and 29 DF, p-value: 1.229e-10Incercam sa mai imbunatatim modelul prin adaugarea cubului lui disp. P-value si multiple r-squared sunt mai bune.
car.order3 <- lm(mpg~disp+I(disp^2)+I(disp^3),data=mtcars)
summary(car.order3)##
## Call:
## lm(formula = mpg ~ disp + I(disp^2) + I(disp^3), data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0896 -1.5653 -0.3619 1.4368 4.7617
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.070e+01 3.809e+00 13.310 1.25e-13 ***
## disp -3.372e-01 5.526e-02 -6.102 1.39e-06 ***
## I(disp^2) 1.109e-03 2.265e-04 4.897 3.68e-05 ***
## I(disp^3) -1.217e-06 2.776e-07 -4.382 0.00015 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.224 on 28 degrees of freedom
## Multiple R-squared: 0.8771, Adjusted R-squared: 0.8639
## F-statistic: 66.58 on 3 and 28 DF, p-value: 7.347e-13b) Transformare logaritmica
In situatiile de modelare statistica, cand avem observatii numerice pozitive, este comun sa aplicam o transformare log pe date pentru a reduce dramatic range-ul datelor si pentru a aduce observatiile extreme mai aproape de centru.
Plot pentru valorile log ale intregilor de la 1 la 1000 (si cele negative) si valorile intregilor. Se poate observa cum valorile log cresc si apoi se aplatizeaza pe masura ce creste valoarea intregului.
plot(1:1000,log(1:1000),type="l",xlab="x",ylab="",ylim=c(-8,8))
lines(1:1000,-log(1:1000),lty=2)
legend("topleft",legend=c("log(x)","-log(x)"),lty=c(1,2))
Ne intoarcem la mtcars si consideram kilometrajul ca functie de cai putere si tipul transmisiei (hp sau am). Construim scatterplot-ul si observam ca punctele sugereaza ca o curba ar fi mai potrivita decat o relatie liniara.
plot(mtcars$hp,mtcars$mpg,pch=19,col=c("black","gray")[factor(mtcars$am)],
xlab="Horsepower",ylab="MPG")
legend("topright",legend=c("auto","man"),col=c("black","gray"),pch=19)
Pentru aceasta folosim transformarea log a cailor putere si tipul transmisiei. Rezultatul confirma efectul semnificativ al cailor putere si al tipului transmisiei asupra kilometrajului. Daca se mentine transmisia constanta, media MPG scade cu 9.24 pentru fiecare unitate a log horsepower. Daca avem transmisie manuala, media MPG creste cu 4.2. Coeficientul de determinare arata ca 82.7% din variatia raspunsului este explicata de aceasta regresie, deci e un model potrivit.
car.log <- lm(mpg~log(hp)+am,data=mtcars)
summary(car.log)##
## Call:
## lm(formula = mpg ~ log(hp) + am, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.9084 -1.7692 -0.1432 1.4032 6.3865
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 63.4842 5.2697 12.047 8.24e-13 ***
## log(hp) -9.2383 1.0439 -8.850 9.78e-10 ***
## am 4.2025 0.9942 4.227 0.000215 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.592 on 29 degrees of freedom
## Multiple R-squared: 0.827, Adjusted R-squared: 0.8151
## F-statistic: 69.31 on 2 and 29 DF, p-value: 8.949e-12plot(mtcars$hp,mtcars$mpg,pch=19,col=c("black","gray")[factor(mtcars$am)],
xlab="Horsepower",ylab="MPG")
legend("topright",legend=c("auto","man"),col=c("black","gray"),pch=19)
hp.seq <- seq(min(mtcars$hp)-20,max(mtcars$hp)+20,length=30)
n <- length(hp.seq)
car.log.pred <- predict(car.log,newdata=data.frame(hp=rep(hp.seq,2),
am=rep(c(0,1),each=n)))
lines(hp.seq,car.log.pred[1:n])
lines(hp.seq,car.log.pred[(n+1):(2*n)],col="gray")
2.10. Exercitii
Exercitiul 21.1
a.
#library("MASS")
?cats
#a
plot(cats$Bwt,cats$Hwt,col=cats$Sex,xlab="Body weight (kg)",ylab="Heart weight (g)")
legend("topleft",legend=c("female","male"),col=c(1,2),pch=c(1,1))
b.
cats.fit <- lm(Hwt~Bwt+Sex,data=cats)
summary(cats.fit)##
## Call:
## lm(formula = Hwt ~ Bwt + Sex, data = cats)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5833 -0.9700 -0.0948 1.0432 5.1016
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.4149 0.7273 -0.571 0.569
## Bwt 4.0758 0.2948 13.826 <2e-16 ***
## SexM -0.0821 0.3040 -0.270 0.788
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.457 on 141 degrees of freedom
## Multiple R-squared: 0.6468, Adjusted R-squared: 0.6418
## F-statistic: 129.1 on 2 and 141 DF, p-value: < 2.2e-16b.i “Mean heart weight” = -0.415 + 4.076* “Body weight” - 0.082* “is male”
Pentru pisicile de acelasi sex, efectul unui kg in plus la greutate inseamna, in medie 4.076 g in plus la greutatea inimii. Pentru pisicile cu aceeasi greutate, greutatea inimii unui mascul este, in medie, cu 0.082 g mai mica decat cea a unei femele. Modelul arata ca efectul greutatii corpului este semnificativ din punct de vedere statistic. Exista dovezi care sugereaza ca greutatea corpului afecteaza greutatea inimii. Totusi, sexul nu este semnificativ, ceea ce inseamna ca includerea variabilei “sex” ca predictor nu este necesara din punct de vedere statistic.
b.ii
names(summary(cats.fit))## [1] "call" "terms" "residuals" "coefficients"
## [5] "aliased" "sigma" "df" "r.squared"
## [9] "adj.r.squared" "fstatistic" "cov.unscaled"summary(cats.fit)$r.squared## [1] 0.6468035Coeficientul de determinare arata ca pentru acest model, 64.5% din variatia greutatii inimii poate fi cuprinsa de regresie.
summary(cats.fit)$fstatistic## value numdf dendf
## 129.1056 2.0000 141.00001-pf(129.1056,df1=2,df2=141)## [1] 0Rezultatul testului omnibus F-test este o p-value mica, 0, ceea ce infirma null hypothesis.
c.
predict(cats.fit,newdata=data.frame(Bwt=3.4,Sex="F"),interval="prediction",level=0.95)## fit lwr upr
## 1 13.44266 10.46904 16.41628d.
plot(cats$Bwt,cats$Hwt,col=cats$Sex,xlab="Body weight (kg)",ylab="Heart weight (g)")
legend("topleft",legend=c("female","male"),col=c(1,2),pch=c(1,1))
Bwt.seq <- seq(min(cats$Bwt)-0.5,max(cats$Bwt)+0.5,length=30)
n <- length(Bwt.seq)
cats.pred <- predict(cats.fit,newdata=data.frame(Bwt=rep(Bwt.seq,2),Sex=rep(c("M","F"),each=n)))
lines(Bwt.seq,cats.pred[1:n],col=2)
lines(Bwt.seq,cats.pred[(n+1):(2*n)])
Cele doua linii suprapuse au panta pozitiva datorita coeficientului “Bwt”, dar foarte aproape una de alta, aratand impactul minimal al variabilei “Sex”.
e.
library("boot")
?nuclear
pairs(nuclear)
f.
nuc.fit1 <- lm(cost~t1+t2,data=nuclear)
summary(nuc.fit1)##
## Call:
## lm(formula = cost ~ t1 + t2, data = nuclear)
##
## Residuals:
## Min 1Q Median 3Q Max
## -273.17 -73.42 -13.40 69.31 360.61
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -242.146 268.020 -0.903 0.37373
## t1 29.908 9.086 3.292 0.00262 **
## t2 4.689 2.945 1.592 0.12224
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 150.1 on 29 degrees of freedom
## Multiple R-squared: 0.272, Adjusted R-squared: 0.2218
## F-statistic: 5.418 on 2 and 29 DF, p-value: 0.01001g.
nuc.fit2 <- lm(cost~t1+t2+date,data=nuclear)
summary(nuc.fit2)##
## Call:
## lm(formula = cost ~ t1 + t2 + date, data = nuclear)
##
## Residuals:
## Min 1Q Median 3Q Max
## -208.63 -90.74 -12.07 59.78 324.19
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9232.833 2974.432 -3.104 0.00434 **
## t1 -5.918 14.281 -0.414 0.68176
## t2 4.639 2.601 1.784 0.08535 .
## date 138.324 45.617 3.032 0.00519 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 132.5 on 28 degrees of freedom
## Multiple R-squared: 0.452, Adjusted R-squared: 0.3933
## F-statistic: 7.698 on 3 and 28 DF, p-value: 0.000667Dupa ce am inclus “date” in modelul liniar s-a redus complet semnificatia statistica a variabilei “t1” din modelul anterior (si-a schimbat semnul). Acest lucru implica faptul ca relatia anterioara pozitiva dintre “t1” si “cost” este explicata mai mult de “date”, care ar trebui folosit intr-un model. Predictorul “t2” ramane nesemnificativ, desi are o p-value mica in acest model.
h.
nuc.fit3 <- lm(cost~date+cap+ne,data=nuclear)
summary(nuc.fit3) ##
## Call:
## lm(formula = cost ~ date + cap + ne, data = nuclear)
##
## Residuals:
## Min 1Q Median 3Q Max
## -154.966 -68.202 -3.614 45.919 285.014
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.458e+03 1.216e+03 -5.310 1.19e-05 ***
## date 9.544e+01 1.773e+01 5.382 9.77e-06 ***
## cap 4.157e-01 9.463e-02 4.393 0.000145 ***
## ne 1.261e+02 4.092e+01 3.083 0.004575 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 99.74 on 28 degrees of freedom
## Multiple R-squared: 0.6895, Adjusted R-squared: 0.6562
## F-statistic: 20.73 on 3 and 28 DF, p-value: 2.827e-07# "cost" = -6458 + 95.4*"date of permit issue" + 0.42*"capacity" + 126.1*"constructed in north-east"
confint(nuc.fit3)## 2.5 % 97.5 %
## (Intercept) -8949.8032112 -3966.9745900
## date 59.1134640 131.7636535
## cap 0.2218791 0.6095524
## ne 42.3145363 209.9430014i.
detroit <- data.frame(Murder=c(8.6,8.9,8.52,8.89,13.07,14.57,21.36,28.03,31.49,37.39,46.26,47.24,52.33),Police=c(260.35,269.8,272.04,272.96,272.51,261.34,268.89,295.99,319.87,341.43,356.59,376.69,390.19),Unemploy=c(11,7,5.2,4.3,3.5,3.2,4.1,3.9,3.6,7.1,8.4,7.7,6.3),Guns=c(178.15,156.41,198.02,222.1,301.92,391.22,665.56,1131.21,837.6,794.9,817.74,583.17,709.59))
detroit## Murder Police Unemploy Guns
## 1 8.60 260.35 11.0 178.15
## 2 8.90 269.80 7.0 156.41
## 3 8.52 272.04 5.2 198.02
## 4 8.89 272.96 4.3 222.10
## 5 13.07 272.51 3.5 301.92
## 6 14.57 261.34 3.2 391.22
## 7 21.36 268.89 4.1 665.56
## 8 28.03 295.99 3.9 1131.21
## 9 31.49 319.87 3.6 837.60
## 10 37.39 341.43 7.1 794.90
## 11 46.26 356.59 8.4 817.74
## 12 47.24 376.69 7.7 583.17
## 13 52.33 390.19 6.3 709.59pairs(detroit)
Numarul politiei pare cea mai sugestiva variabila pentru a prezice numarul de crime.
j.
murd.fit <- lm(Murder~Police+Unemploy+Guns,data=detroit)
summary(murd.fit)##
## Call:
## lm(formula = Murder ~ Police + Unemploy + Guns, data = detroit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.8422 -1.9451 0.2012 0.9172 4.6694
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -68.852509 5.862631 -11.744 9.25e-07 ***
## Police 0.280799 0.024657 11.388 1.20e-06 ***
## Unemploy 0.147248 0.408235 0.361 0.72665
## Guns 0.014177 0.003538 4.007 0.00308 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.89 on 9 degrees of freedom
## Multiple R-squared: 0.9767, Adjusted R-squared: 0.9689
## F-statistic: 125.6 on 3 and 9 DF, p-value: 1.158e-07summary(murd.fit)$r.squared## [1] 0.9766753“Mean murders” = -68.85 + 0.281* “no. of police” + 0.147* “unemployment” + 0.014*“no. of gun licenses”.
Pentru “no. of gun licenses” si “unemployment”, fiecare politie in plus per 10000 populatie inseamna o crestere medie de 0.28 crime per 10000 populatie. Pentru “no. of police” si “unemployment”, fiecare licenta de arma in plus per 10000 populatie inseamna o crestere medie de 0.014 crime per 10000 populatie. Pentru “no. of gun licenses” si “no. of police”, fiecare procent de somaj inseamna o crestere medie de 0.147 crime per 10000 populatie. Nu putem justifica ca vreuna din relatii e cauzala, bazandu-ne doar pe setul de date si analiza. Este greu de justificat ca avand mai multi membri in politie creste numarul de crime, de exemplu.
k.
summary(murd.fit)$r.squared## [1] 0.9766753Aprox. 97.67% din variabilitatea nr. mediu de crime este explicata de modelul cu 3 predictori.
murd.fit2 <- lm(Murder~Police+Guns,data=detroit)
summary(murd.fit2)##
## Call:
## lm(formula = Murder ~ Police + Guns, data = detroit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.9424 -2.1068 0.2775 0.9614 4.6408
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -69.002919 5.587647 -12.349 2.23e-07 ***
## Police 0.285048 0.020697 13.772 7.92e-08 ***
## Guns 0.013636 0.003062 4.453 0.00123 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.761 on 10 degrees of freedom
## Multiple R-squared: 0.9763, Adjusted R-squared: 0.9716
## F-statistic: 206.3 on 2 and 10 DF, p-value: 7.417e-09summary(murd.fit2)$r.squared## [1] 0.9763381Coeficientul de determinare nu se schimba foarte mult, fiind 97.63% pentru modelul cu 2 predictori.
l.
predict(murd.fit2,newdata=data.frame(Police=c(300,300),Guns=c(500,0)),interval="confidence",level=0.99)## fit lwr upr
## 1 23.32948 20.88251 25.77645
## 2 16.51159 10.90530 22.11787Exercitiul 21.2
a.
gal <- data.frame(d=c(573,534,495,451,395,337,253),h=c(1,0.8,0.6,0.45,0.3,0.2,0.1))
plot(gal$d~gal$h,pch=19,xlab="Initial height",ylab="Distance traveled")
b.i
gal.fit.order2 <- lm(d~h+I(h^2),data=gal)
summary(gal.fit.order2)##
## Call:
## lm(formula = d ~ h + I(h^2), data = gal)
##
## Residuals:
## 1 2 3 4 5 6 7
## 8.458 -12.607 -6.177 1.940 13.523 9.170 -14.308
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 199.91 16.76 11.928 0.000283 ***
## h 708.32 74.82 9.467 0.000695 ***
## I(h^2) -343.69 66.78 -5.147 0.006760 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.64 on 4 degrees of freedom
## Multiple R-squared: 0.9903, Adjusted R-squared: 0.9855
## F-statistic: 205 on 2 and 4 DF, p-value: 9.333e-05b.ii
gal.fit.order3 <- lm(d~h+I(h^2)+I(h^3),data=gal)
summary(gal.fit.order3)##
## Call:
## lm(formula = d ~ h + I(h^2) + I(h^3), data = gal)
##
## Residuals:
## 1 2 3 4 5 6 7
## -0.84138 2.32159 -0.08044 -4.46885 1.89175 3.58091 -2.40359
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 155.776 8.326 18.710 0.000333 ***
## h 1115.298 65.671 16.983 0.000445 ***
## I(h^2) -1244.943 138.425 -8.994 0.002902 **
## I(h^3) 547.710 83.273 6.577 0.007150 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.011 on 3 degrees of freedom
## Multiple R-squared: 0.9994, Adjusted R-squared: 0.9987
## F-statistic: 1595 on 3 and 3 DF, p-value: 2.662e-05gal.fit.order4 <- lm(d~h+I(h^2)+I(h^3)+I(h^4),data=gal)
summary(gal.fit.order4)##
## Call:
## lm(formula = d ~ h + I(h^2) + I(h^3) + I(h^4), data = gal)
##
## Residuals:
## 1 2 3 4 5 6 7
## 0.1708 -0.9279 2.3183 -2.3092 0.2576 0.9338 -0.4433
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 138.295 9.066 15.254 0.00427 **
## h 1346.071 106.071 12.690 0.00615 **
## I(h^2) -2116.913 379.271 -5.582 0.03063 *
## I(h^3) 1766.391 518.566 3.406 0.07644 .
## I(h^4) -561.014 237.498 -2.362 0.14201
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.523 on 2 degrees of freedom
## Multiple R-squared: 0.9998, Adjusted R-squared: 0.9995
## F-statistic: 3024 on 4 and 2 DF, p-value: 0.0003306Modelul de ordin 3 este semnificativ, cu un coeficient de determinare mai bun, spre deosebire de modelul de ordin 4.
c. Din cele 3 modele, functie cubica pare cea mai buna, relatia dintre “distance traveled” si “initial height” parand cubica. Modelul de ordin 2 pare prea simplu, in timp ce modelul de ordin 4 pare prea complex.
plot(gal$d~gal$h,pch=19,xlab="Initial height",ylab="Distance traveled")
hseq <- seq(0.05,1.05,length=30)
gal.pred <- predict(gal.fit.order3,newdata=data.frame(h=hseq),interval="confidence",level=0.9)
lines(hseq,gal.pred[,1])
lines(hseq,gal.pred[,2],lty=2)
lines(hseq,gal.pred[,3],lty=2)
d.
library("faraway")##
## Attaching package: 'faraway'## The following objects are masked from 'package:boot':
##
## logit, melanoma?trees
plot(trees$Volume~trees$Girth,pch=19,xlab="Girth",ylab="Volume")
e.
tree.fit1 <- lm(Volume~Girth+I(Girth^2),trees)
summary(tree.fit1) ## "Mean volume" = 10.79 - 2.09*"girth" + 0.254*"girth^2"##
## Call:
## lm(formula = Volume ~ Girth + I(Girth^2), data = trees)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.4889 -2.4293 -0.3718 2.0764 7.6447
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.78627 11.22282 0.961 0.344728
## Girth -2.09214 1.64734 -1.270 0.214534
## I(Girth^2) 0.25454 0.05817 4.376 0.000152 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.335 on 28 degrees of freedom
## Multiple R-squared: 0.9616, Adjusted R-squared: 0.9588
## F-statistic: 350.5 on 2 and 28 DF, p-value: < 2.2e-16tree.fit2 <- lm(log(Volume)~log(Girth),trees)
summary(tree.fit2) ## "Mean log(volume)" = -2.35 + 2.20*"log(girth)"##
## Call:
## lm(formula = log(Volume) ~ log(Girth), data = trees)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.205999 -0.068702 0.001011 0.072585 0.247963
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.35332 0.23066 -10.20 4.18e-11 ***
## log(Girth) 2.19997 0.08983 24.49 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.115 on 29 degrees of freedom
## Multiple R-squared: 0.9539, Adjusted R-squared: 0.9523
## F-statistic: 599.7 on 1 and 29 DF, p-value: < 2.2e-16Coeficientii de determinare sunt similari, modelul de ordin 4 fiind putin mai mare. Ambele indica un efect pozitiv semnificativ din punct de vedere statistic asupra girth on volume.
f.
gseq <- seq(7,21,length=30)
tree.pred1 <- predict(tree.fit1,newdata=data.frame(Girth=gseq),interval="prediction")
tree.pred2 <- predict(tree.fit2,newdata=data.frame(Girth=gseq),interval="prediction")
plot(trees$Volume~trees$Girth,pch=19,xlab="Girth",ylab="Volume")
lines(gseq,tree.pred1[,1],lwd=2)
lines(gseq,tree.pred1[,2])
lines(gseq,tree.pred1[,3])
lines(gseq,exp(tree.pred2[,1]),lwd=2,lty=2)
lines(gseq,exp(tree.pred2[,2]),lty=2)
lines(gseq,exp(tree.pred2[,3]),lty=2)
legend("topleft",legend=c("quadratic","logged"),lty=1:2)
Valorile fitted ale modelelor sunt foarte similare. Totusi, modelul de ordin 4 are limite mult mai mari pentru valori girth mici decat pentru valori girth mari. De cealalta parte, limitele pentru modelul cu log sunt substantial mai mari decat acelea ale modelului patratic la valori girth mai mari.
g.
#library("MASS")
car.fit <- lm(mpg~wt+hp+disp,data=mtcars)
summary(car.fit)##
## Call:
## lm(formula = mpg ~ wt + hp + disp, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.891 -1.640 -0.172 1.061 5.861
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.105505 2.110815 17.579 < 2e-16 ***
## wt -3.800891 1.066191 -3.565 0.00133 **
## hp -0.031157 0.011436 -2.724 0.01097 *
## disp -0.000937 0.010350 -0.091 0.92851
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.639 on 28 degrees of freedom
## Multiple R-squared: 0.8268, Adjusted R-squared: 0.8083
## F-statistic: 44.57 on 3 and 28 DF, p-value: 8.65e-11h.
car.fit <- lm(I(1/mpg)~wt+hp+disp,data=mtcars)
summary(car.fit)##
## Call:
## lm(formula = I(1/mpg) ~ wt + hp + disp, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0163719 -0.0043511 0.0008672 0.0032544 0.0133345
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.496e-03 5.322e-03 1.784 0.08521 .
## wt 9.469e-03 2.688e-03 3.522 0.00149 **
## hp 5.864e-05 2.883e-05 2.034 0.05155 .
## disp 2.456e-05 2.609e-05 0.941 0.35472
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.006653 on 28 degrees of freedom
## Multiple R-squared: 0.8518, Adjusted R-squared: 0.8359
## F-statistic: 53.63 on 3 and 28 DF, p-value: 9.94e-12Ambele modele au nivele similare de semnificatie pentru cei 3 predictori. Totusi, se observa o imbunatatire a coeficientului de determinare pentru ultimul model bazat pe variabila de raspuns GPM = 1/MPG.
Exercitiul 21.3
a.
cat.fit <- lm(Hwt~Bwt*Sex,data=cats)
summary(cat.fit)##
## Call:
## lm(formula = Hwt ~ Bwt * Sex, data = cats)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.7728 -1.0118 -0.1196 0.9272 4.8646
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.9813 1.8428 1.618 0.107960
## Bwt 2.6364 0.7759 3.398 0.000885 ***
## SexM -4.1654 2.0618 -2.020 0.045258 *
## Bwt:SexM 1.6763 0.8373 2.002 0.047225 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.442 on 140 degrees of freedom
## Multiple R-squared: 0.6566, Adjusted R-squared: 0.6493
## F-statistic: 89.24 on 3 and 140 DF, p-value: < 2.2e-16Modelul main-effects-only avea un efect usor negativ al variabilei “sex male” si nu era semnificativ. In aceasta versiune, efectul “sex male” e mai mare si p-value este mult mai mica, de unde rezulta dovada slaba pentru semnificatie.
b.
plot(cats$Bwt,
cats$Hwt,
col=cats$Sex,
ylab="Heart weight (g)",
xlab="Body weight (kg)")
legend("topleft",
legend=c("Female","Male"),
col=1:2,pch=1)
cat.coefs <- coef(cat.fit)
abline(coef=cat.coefs[1:2])
abline(coef=c(sum(cat.coefs[c(1,3)]),sum(cat.coefs[c(2,4)])),col=2)
Liniile modelului nu mai sunt paralele.
c.
predict(cat.fit,
newdata=data.frame(Bwt=3.4,
Sex="F"),
interval="prediction",
level=0.95)## fit lwr upr
## 1 11.94512 8.651786 15.23845Greutatea inimii lui Sigma prezisa de noul model este cu aproximativ 1.5g mai mica decat cea prezisa de modelul anterior main-effects-only.Intervalul de predictie este, de asemenea, mai jos, dar este mai larg decat cel anterior.
d.
#library("faraway")
tree.fit1 <- lm(Volume~Girth+Height,data=trees)
summary(tree.fit1)##
## Call:
## lm(formula = Volume ~ Girth + Height, data = trees)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.4065 -2.6493 -0.2876 2.2003 8.4847
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -57.9877 8.6382 -6.713 2.75e-07 ***
## Girth 4.7082 0.2643 17.816 < 2e-16 ***
## Height 0.3393 0.1302 2.607 0.0145 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.882 on 28 degrees of freedom
## Multiple R-squared: 0.948, Adjusted R-squared: 0.9442
## F-statistic: 255 on 2 and 28 DF, p-value: < 2.2e-16tree.fit2 <- lm(Volume~Girth*Height,data=trees)
summary(tree.fit2)##
## Call:
## lm(formula = Volume ~ Girth * Height, data = trees)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.5821 -1.0673 0.3026 1.5641 4.6649
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 69.39632 23.83575 2.911 0.00713 **
## Girth -5.85585 1.92134 -3.048 0.00511 **
## Height -1.29708 0.30984 -4.186 0.00027 ***
## Girth:Height 0.13465 0.02438 5.524 7.48e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.709 on 27 degrees of freedom
## Multiple R-squared: 0.9756, Adjusted R-squared: 0.9728
## F-statistic: 359.3 on 3 and 27 DF, p-value: < 2.2e-16e.
tree.fit3 <- lm(log(Volume)~log(Girth)+log(Height),data=trees)
summary(tree.fit3)##
## Call:
## lm(formula = log(Volume) ~ log(Girth) + log(Height), data = trees)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.168561 -0.048488 0.002431 0.063637 0.129223
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.63162 0.79979 -8.292 5.06e-09 ***
## log(Girth) 1.98265 0.07501 26.432 < 2e-16 ***
## log(Height) 1.11712 0.20444 5.464 7.81e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08139 on 28 degrees of freedom
## Multiple R-squared: 0.9777, Adjusted R-squared: 0.9761
## F-statistic: 613.2 on 2 and 28 DF, p-value: < 2.2e-16tree.fit4 <- lm(log(Volume)~log(Girth)*log(Height),data=trees)
summary(tree.fit4)##
## Call:
## lm(formula = log(Volume) ~ log(Girth) * log(Height), data = trees)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.165941 -0.048613 0.006384 0.062204 0.132295
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.6869 7.6996 -0.479 0.636
## log(Girth) 0.7942 3.0910 0.257 0.799
## log(Height) 0.4377 1.7788 0.246 0.808
## log(Girth):log(Height) 0.2740 0.7124 0.385 0.704
##
## Residual standard error: 0.08265 on 27 degrees of freedom
## Multiple R-squared: 0.9778, Adjusted R-squared: 0.9753
## F-statistic: 396.4 on 3 and 27 DF, p-value: < 2.2e-16Efectul este foarte semnificativ in modelul netransformat de la punctul d, dar complet nesemnificativ dupa transformarile logaritmice a tuturor variabilelor. Acest lucru sugereaza ca relatiile liniare nu sunt cel mai bun mod de a modela aceste date si ca trebuie sa luam in calcul curbura suprafetei de raspuns fie prin transformarea datelor, fie prin includerea unei interactiuni bidirectionale intre cei doi predictori netransformati. Este dificil de ales una dintre aceste solutii, deoarece trebuie sa stim mai multe atat despre natura datelor in context, cat si despre scopul modelului.
f.
car.fit <- lm(mpg~factor(cyl)*hp+wt,data=mtcars)
summary(car.fit)##
## Call:
## lm(formula = mpg ~ factor(cyl) * hp + wt, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.1864 -1.4098 -0.4022 1.0186 4.3920
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 41.87732 3.23293 12.953 1.37e-12 ***
## factor(cyl)6 -9.98213 5.76950 -1.730 0.095931 .
## factor(cyl)8 -11.72793 4.22507 -2.776 0.010276 *
## hp -0.09947 0.03487 -2.853 0.008576 **
## wt -3.05994 0.68275 -4.482 0.000143 ***
## factor(cyl)6:hp 0.07809 0.05236 1.492 0.148335
## factor(cyl)8:hp 0.08602 0.03703 2.323 0.028601 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.3 on 25 degrees of freedom
## Multiple R-squared: 0.8826, Adjusted R-squared: 0.8544
## F-statistic: 31.32 on 6 and 25 DF, p-value: 1.831e-10g.
coef(car.fit)## (Intercept) factor(cyl)6 factor(cyl)8 hp wt
## 41.87732085 -9.98213264 -11.72792959 -0.09946598 -3.05993524
## factor(cyl)6:hp factor(cyl)8:hp
## 0.07808919 0.08602496Efectul este intre un predictor continuu (hp) si unul categoric (factor(cyl)). Oricare dintre cei doi coeficienti estimati poate fi interpretat ca schimbarea pantei a lui hp pentru fiecare dintre nivelurile non-reference a variabilei factor(cyl).
coef(car.fit)[4] # cand o masina are 4 cilindri (nivel de referinta), panta pentru hp este ## hp
## -0.09946598# -0.0995
coef(car.fit)[4] + coef(car.fit)[6] # cand masina are 6 cilindri, panta pentru hp este ## hp
## -0.02137679# -0.0995 + 0.0781 = -0.0214
coef(car.fit)[4] + coef(car.fit)[7] # cand masina are 8 cilindri, panta pentru hp este## hp
## -0.01344103# -0.0995 + 0.0860 = -0.0135Acest model sugereaza ca cu cat hp creste, cu atat media MPG scade (pentru wt fixat). Totusi, in comparatie cu masinile cu 4 cilindri, se estimeaza ca media MPG sa scada intr-un ritm mai lent cu valoare hp care creste pentru masinile cu 6 si 8 cilindri.
h.i
predict(car.fit,
newdata=data.frame(wt=c(2.1,3.9,2.9),
hp=c(100,210,200),
cyl=c(4,8,6)),
interval="confidence",
level=0.95)## fit lwr upr
## 1 25.50486 23.57668 27.43304
## 2 15.39303 14.11928 16.66678
## 3 18.74602 12.29560 25.19644Prima masina este singura care are o estimare a mediei MPG mai mare decat cererea de 25, deci aceasta ar fi alegerea initiala.
h.ii
Desi estimarea pentru masina 3 este mult mai mica decat 25, uitandu-ne la intervalele de incredere putem vedea ca intervalul pentru masina 3 include 25. Deci, putem spune ca suntem siguri in proportie de 95% ca medie MPG adevarata a unei masini precum masina 3 se afla undeva in acel interval. Modelul nu sugereaza nicio dovada impotriva ipotezei conform careia media MPG adevarata a unei astfel de masini este egala cu 25.
3. Selectarea și cercetarea modelelor de regresii liniare
3.1 Complexitate vs. Acuratețe
Concepte esențiale în validarea modelelor de regresie liniară:
- selectarea unui model potrivit analizei
- validarea ipotezelor de lucru
În vederea selectării unui model, trebuie să se ajungă la un compromis între complexitate - direct proporțională cu numărul de variabile luate în considerare - și perfecționarea modelului. În acest sens este utilizat termenul de Principiu al parciomoniei ce presupune obținerea unui model de complexitate cât mai scăzută, dar fără a sacrifica prea mult din precizia modelului.
3.2 Algoritmi de selecție a modelelor
Exista multe metode de selecție a modelelor și nicio abordare nu este universal valabilă pentru fiecare model de regresie. Algoritmii de selecție diferiți pot avea ca rezultat diferite modele finale. În multe cazuri, cercetătorii vor avea suplimentar informații sau cunoștințe despre problema care influențează decizia — de exemplu, că anumiți predictori trebuie să fie întotdeauna incluși sau că nu are sens să le includă.
3.2.1 Nested Comparisons: The Partial F-Test
Testul-F Parțial este o metodă de comparare a două modele, dintre care modelul mai puțin complex este o versiune redusă a celuilalt.
\[ \hat{y}_{redu} = \hat{\beta_0} + \hat{\beta_{1}}\ x_1 + \hat{\beta_{2}}\ x_2 + ... + \hat{\beta_{p}}\ x_p \] \[ \hat{y}_{full} = \hat{\beta_0} + \hat{\beta_{1}}\ x_1 + \hat{\beta_{2}}\ x_2 + ... + \hat{\beta_{p}}\ x_p + ... + \hat{\beta_{q}}\ x_q \]
Modelul \(\hat{y}_{redu}\) are \(p\) termeni plus cel de interceptare. Modelul \(\hat{y}_{full}\), care în include pe cel \(\hat{y}_{redu}\) are \(q\) - \(p\) termeni în plus. Problema este pusă în felul următor: este noua îmbunătățire a primului model suficient de bună?
Testul-F Parțial se calculează cu următoarea formulă:
- \(n\) - numărul de variabile ale eșantionului de date
- \(p\) - numărul de variabile al modelului redus
- \(q\) - numărul de variabile al modelului complex
Exemplu: Pentru compararea a două modele de regresii liniare diferite, se folosește funcția \(anova\) ce admite ca parametrii regresiile care vor fi comparate. Am folosit setul survey ce conține datele a 237 de studenți cu privire la gen, lungimea palmelor, înalțimea, frecvența fumatului etc. Testăm dacă adăugarea factorului fumat la un model de regresie deja generat care prezice înălțimea studenului îmbunătățește semnificativ modelul:
#library(MASS)
#data(survey)
survmult <- lm(Height ~ Wr.Hnd + Sex, data = survey) # model care ia în calcul lungimea palmei cu care scrie studentul și genul
survmult2 <- lm(Height ~ Wr.Hnd + Sex + Smoke, data = survey) # model care ia în calcul lungimea palmei cu care scrie studentul, genul și factorul fumat
anova(survmult, survmult2) # Se compară cele 2 modele## Analysis of Variance Table
##
## Model 1: Height ~ Wr.Hnd + Sex
## Model 2: Height ~ Wr.Hnd + Sex + Smoke
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 204 9959.2
## 2 201 9914.3 3 44.876 0.3033 0.823Cum \(p-value\) = \(0.823\) >>>> \(0.05\) (valoare de prag aleasă), acest lucru indică faptul că adăgarea factorului fumat la forma redusă a modelului nu oferă îmbunătățiri suficient de bune pentru prezicerea înălțimii.
3.2.2 Forward Selection
Testul-F Parțial reprezintă o modalitate naturală de cercetare, dar acest lucru poate deveni dificil în momentul în care există o multitudine de combinații de variabile pentru un model. Astfel a apărut Forward-Selection care pornește de la un model trivial, după care se execută o serie de teste pentru a vedea care variabilă îmbunătățește cel mai mult modelul. Acesta se actualizează și se repetă procesul până când nu există termeni care să perfecționeze semnificativ modelul.
Exemplu: Pentru a ilustra funcționarea algoritmului, se va folosi setul de date nuclear ce conține date despre centrale cu reactoare cu apă ușoară în vederea prezicerii costului de construcție.
Se folosesc funcțiile:
- \(add1\) pentru a calcula toți termenii unici din argumentul scope care pot fi adăugați la model, îi antrenează și generează un tabel cu noile rezultate, având următorii parametri:
- object - obiectul cu regresia liniară deja construită
- scope - formula cu termenii care pot fi adăugați la regresia liniară
- test - tipul de test efectuat cu modelul original
- \(update\) ce va actualiza și reantrena un model, având parametrii:
- object - modelul de regresie liniară care va fi actualizat
- formula - noua formulă cu care va fi actualizat modelul
#library(boot)
nuc.0 <- lm(cost ~ 1, data = nuclear)
summary(nuc.0)##
## Call:
## lm(formula = cost ~ 1, data = nuclear)
##
## Residuals:
## Min 1Q Median 3Q Max
## -254.05 -151.24 -13.46 150.40 419.68
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 461.56 30.07 15.35 4.95e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 170.1 on 31 degrees of freedomadd1(nuc.0, scope = .~. + date + t1 + t2 + cap + pr +ne + ct + bw + cum.n + pt, test = "F")## Single term additions
##
## Model:
## cost ~ 1
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 897172 329.72
## date 1 334335 562837 316.80 17.8205 0.0002071 ***
## t1 1 186984 710189 324.24 7.8986 0.0086296 **
## t2 1 27 897145 331.72 0.0009 0.9760597
## cap 1 199673 697499 323.66 8.5881 0.0064137 **
## pr 1 9037 888136 331.40 0.3052 0.5847053
## ne 1 128641 768531 326.77 5.0216 0.0325885 *
## ct 1 43042 854130 330.15 1.5118 0.2284221
## bw 1 16205 880967 331.14 0.5519 0.4633402
## cum.n 1 67938 829234 329.20 2.4579 0.1274266
## pt 1 305334 591839 318.41 15.4772 0.0004575 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1nuc.1 <- update(nuc.0, formula = .~. + date)
summary(nuc.1)##
## Call:
## lm(formula = cost ~ date, data = nuclear)
##
## Residuals:
## Min 1Q Median 3Q Max
## -176.00 -105.27 -25.24 58.63 359.46
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6553.57 1661.96 -3.943 0.000446 ***
## date 102.29 24.23 4.221 0.000207 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 137 on 30 degrees of freedom
## Multiple R-squared: 0.3727, Adjusted R-squared: 0.3517
## F-statistic: 17.82 on 1 and 30 DF, p-value: 0.0002071add1(nuc.1, scope = .~. + date + t1 + t2 +cap + pr + ne + ct + bw + cum.n + pt, test="F")## Single term additions
##
## Model:
## cost ~ date
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 562837 316.80
## t1 1 15322 547515 317.92 0.8115 0.3750843
## t2 1 68161 494676 314.67 3.9959 0.0550606 .
## cap 1 189732 373105 305.64 14.7471 0.0006163 ***
## pr 1 4027 558810 318.57 0.2090 0.6509638
## ne 1 92256 470581 313.07 5.6854 0.0238671 *
## ct 1 54794 508043 315.52 3.1277 0.0874906 .
## bw 1 1240 561597 318.73 0.0640 0.8020147
## cum.n 1 4658 558179 318.53 0.2420 0.6264574
## pt 1 90587 472250 313.18 5.5628 0.0252997 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1nuc.2 <- update(nuc.1, formula = .~. + cap)
summary(nuc.2)##
## Call:
## lm(formula = cost ~ date + cap, data = nuclear)
##
## Residuals:
## Min 1Q Median 3Q Max
## -146.351 -88.244 -2.202 58.457 267.741
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6790.8792 1377.6683 -4.929 3.09e-05 ***
## date 100.7764 20.0696 5.021 2.39e-05 ***
## cap 0.4132 0.1076 3.840 0.000616 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 113.4 on 29 degrees of freedom
## Multiple R-squared: 0.5841, Adjusted R-squared: 0.5555
## F-statistic: 20.37 on 2 and 29 DF, p-value: 2.984e-06add1(nuc.2, scope = .~. + date + t1 + t2 + cap + pr + ne + ct + bw + cum.n + pt, test="F")## Single term additions
##
## Model:
## cost ~ date + cap
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 373105 305.64
## t1 1 957 372148 307.56 0.0720 0.790387
## t2 1 13356 359750 306.48 1.0395 0.316670
## pr 1 18453 354652 306.02 1.4569 0.237521
## ne 1 94537 278569 298.29 9.5022 0.004575 **
## ct 1 48957 324148 303.14 4.2289 0.049169 *
## bw 1 7505 365601 306.99 0.5748 0.454705
## cum.n 1 28062 345044 305.14 2.2772 0.142493
## pt 1 95917 277189 298.13 9.6889 0.004243 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1nuc.3 <- update(nuc.2, formula = .~. + pt)
summary(nuc.3)##
## Call:
## lm(formula = cost ~ date + cap + pt, data = nuclear)
##
## Residuals:
## Min 1Q Median 3Q Max
## -190.11 -56.21 -13.90 35.57 237.92
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.553e+03 1.406e+03 -3.238 0.003091 **
## date 6.852e+01 2.043e+01 3.355 0.002296 **
## cap 4.191e-01 9.441e-02 4.439 0.000128 ***
## pt -1.628e+02 5.229e+01 -3.113 0.004243 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 99.5 on 28 degrees of freedom
## Multiple R-squared: 0.691, Adjusted R-squared: 0.6579
## F-statistic: 20.88 on 3 and 28 DF, p-value: 2.64e-07add1(nuc.3, scope = .~. + date + t1 + t2 + cap + pr + ne + ct + bw + cum.n + pt, test="F")## Single term additions
##
## Model:
## cost ~ date + cap + pt
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 277189 298.13
## t1 1 12 277177 300.13 0.0011 0.97352
## t2 1 10429 266760 298.91 1.0555 0.31335
## pr 1 6017 271171 299.43 0.5991 0.44564
## ne 1 54572 222617 293.12 6.6187 0.01591 *
## ct 1 18119 259069 297.97 1.8884 0.18068
## bw 1 654 276535 300.06 0.0639 0.80241
## cum.n 1 463 276726 300.08 0.0451 0.83334
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1nuc.4 <- update(nuc.3, formula = .~. + ne)
summary(nuc.4)##
## Call:
## lm(formula = cost ~ date + cap + pt + ne, data = nuclear)
##
## Residuals:
## Min 1Q Median 3Q Max
## -157.894 -38.424 -2.493 35.363 267.445
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.756e+03 1.286e+03 -3.699 0.000975 ***
## date 7.102e+01 1.867e+01 3.804 0.000741 ***
## cap 4.198e-01 8.616e-02 4.873 4.29e-05 ***
## pt -1.289e+02 4.950e+01 -2.605 0.014761 *
## ne 9.940e+01 3.864e+01 2.573 0.015908 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 90.8 on 27 degrees of freedom
## Multiple R-squared: 0.7519, Adjusted R-squared: 0.7151
## F-statistic: 20.45 on 4 and 27 DF, p-value: 7.507e-08add1(nuc.4, scope = .~. + date + t1 + t2 + cap + pr + ne + ct + bw + cum.n + pt, test="F")## Single term additions
##
## Model:
## cost ~ date + cap + pt + ne
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 222617 293.12
## t1 1 107.0 222510 295.10 0.0125 0.9118
## t2 1 19229.9 203387 292.23 2.4583 0.1290
## pr 1 5230.8 217386 294.36 0.6256 0.4361
## ct 1 15764.7 206852 292.77 1.9815 0.1711
## bw 1 448.0 222169 295.06 0.0524 0.8207
## cum.n 1 13819.9 208797 293.07 1.7209 0.2010summary(nuc.4)##
## Call:
## lm(formula = cost ~ date + cap + pt + ne, data = nuclear)
##
## Residuals:
## Min 1Q Median 3Q Max
## -157.894 -38.424 -2.493 35.363 267.445
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.756e+03 1.286e+03 -3.699 0.000975 ***
## date 7.102e+01 1.867e+01 3.804 0.000741 ***
## cap 4.198e-01 8.616e-02 4.873 4.29e-05 ***
## pt -1.289e+02 4.950e+01 -2.605 0.014761 *
## ne 9.940e+01 3.864e+01 2.573 0.015908 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 90.8 on 27 degrees of freedom
## Multiple R-squared: 0.7519, Adjusted R-squared: 0.7151
## F-statistic: 20.45 on 4 and 27 DF, p-value: 7.507e-08La fiecare pas se generează un nou raport cu privire la adăugarea câte unei variabile. Pentru fiecare raport, se alege p-value-ul cel mai mic și se actualizează noul model. Se repetă pașii cât timp obțin o valoare suficient de bună pentru \(p-value\) (aleasă, de obicei, cu o margine superioară de \(0.01 - 0.05\)), altfel acesta se oprește și am obținut modelul final.
De menționat că există o anumită doză de subiectivism pentru această metodă deoarece presupune compararea și alegerea valori care pot fi destul de apropiate. De exemplu, la primul pas, se putea alege variabila \(pt\) în detrimentul lui \(date\) și s-ar fi ajuns la un model diferit în final.
3.2.3 Backward Selection
Backward Selection are un comportament asemănător cu Forward Selection, cu excepția faptului că acesta pornește de la cel mai complex model și elimină succesiv câte o variabilă. Alegerea dintre Forward și Backward Selection se face în funcție de situație. De exemplu, dacă cel mai complex model nu este cunoscut sau este prea greu de obținut, Forward Selection este preferat deoarece pornește de la un model simplu, trivial.
Exemplu: Pentru a ilustra funcționarea algoritmului, se va folosi, în continuare, setul de date nuclear.
Adițional se folosește funcția \(drop1\) care se comportă asemănător cu funcția \(add1\), doar că aceasta calculeaza toți termenii unici care pot fi eliminați din modelul dat ca parametru și generează un tabel cu noile rezultate. Parametrii sunt:
- object - obiectul cu regresia liniară deja construită
- scope - formula cu termenii care pot fi eliminați din regresia liniară
- test - tipul de test efectuat cu modelul original
# se pornește de la modelul complet, conținând toți termenii
nuc.0 <- lm(cost ~ date + t1 + t2 + cap + pr + ne + ct + bw + cum.n + pt, data = nuclear)
summary(nuc.0)##
## Call:
## lm(formula = cost ~ date + t1 + t2 + cap + pr + ne + ct + bw +
## cum.n + pt, data = nuclear)
##
## Residuals:
## Min 1Q Median 3Q Max
## -128.608 -46.736 -2.668 39.782 180.365
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.135e+03 2.788e+03 -2.918 0.008222 **
## date 1.155e+02 4.226e+01 2.733 0.012470 *
## t1 5.928e+00 1.089e+01 0.545 0.591803
## t2 4.571e+00 2.243e+00 2.038 0.054390 .
## cap 4.217e-01 8.844e-02 4.768 0.000104 ***
## pr -8.112e+01 4.077e+01 -1.990 0.059794 .
## ne 1.375e+02 3.869e+01 3.553 0.001883 **
## ct 4.327e+01 3.431e+01 1.261 0.221008
## bw -8.238e+00 5.188e+01 -0.159 0.875354
## cum.n -6.989e+00 3.822e+00 -1.829 0.081698 .
## pt -1.925e+01 6.367e+01 -0.302 0.765401
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 82.83 on 21 degrees of freedom
## Multiple R-squared: 0.8394, Adjusted R-squared: 0.763
## F-statistic: 10.98 on 10 and 21 DF, p-value: 2.844e-06drop1(nuc.0, test = "F")## Single term deletions
##
## Model:
## cost ~ date + t1 + t2 + cap + pr + ne + ct + bw + cum.n + pt
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 144065 291.19
## date 1 51230 195295 298.93 7.4677 0.0124702 *
## t1 1 2034 146099 289.64 0.2965 0.5918028
## t2 1 28481 172546 294.97 4.1517 0.0543902 .
## cap 1 155943 300008 312.67 22.7314 0.0001039 ***
## pr 1 27161 171226 294.72 3.9592 0.0597943 .
## ne 1 86581 230646 304.25 12.6207 0.0018835 **
## ct 1 10915 154980 291.53 1.5911 0.2210075
## bw 1 173 144238 289.23 0.0252 0.8753538
## cum.n 1 22939 167004 293.92 3.3438 0.0816977 .
## pt 1 627 144692 289.33 0.0914 0.7654015
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1nuc.1 <- update(nuc.0, .~. - bw)
drop1(nuc.1, test = "F")## Single term deletions
##
## Model:
## cost ~ date + t1 + t2 + cap + pr + ne + ct + cum.n + pt
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 144238 289.23
## date 1 55942 200180 297.72 8.5326 0.007913 **
## t1 1 3124 147362 287.92 0.4765 0.497245
## t2 1 30717 174955 293.41 4.6852 0.041546 *
## cap 1 159976 304214 311.11 24.4005 6.098e-05 ***
## pr 1 27140 171377 292.75 4.1395 0.054122 .
## ne 1 86408 230646 302.25 13.1795 0.001479 **
## ct 1 11815 156053 289.75 1.8021 0.193153
## cum.n 1 24048 168286 292.17 3.6680 0.068557 .
## pt 1 930 145168 287.44 0.1419 0.710039
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1nuc.2 <- update(nuc.1, .~. - pt)
drop1(nuc.2, test = "F")## Single term deletions
##
## Model:
## cost ~ date + t1 + t2 + cap + pr + ne + ct + cum.n
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 145168 287.44
## date 1 93677 238845 301.37 14.8419 0.0008108 ***
## t1 1 2853 148021 286.06 0.4521 0.5080427
## t2 1 31312 176480 291.69 4.9610 0.0359961 *
## cap 1 165594 310762 309.79 26.2363 3.447e-05 ***
## pr 1 29984 175152 291.45 4.7506 0.0397817 *
## ne 1 113659 258827 303.94 18.0077 0.0003069 ***
## ct 1 15602 160770 288.70 2.4719 0.1295556
## cum.n 1 48482 193650 294.66 7.6814 0.0108561 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1nuc.3 <- update(nuc.2, .~. - t1)
drop1(nuc.3, test = "F")## Single term deletions
##
## Model:
## cost ~ date + t2 + cap + pr + ne + ct + cum.n
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 148021 286.06
## date 1 369019 517040 324.08 59.8322 5.697e-08 ***
## t2 1 28539 176560 289.70 4.6272 0.0417479 *
## cap 1 162846 310867 307.80 26.4036 2.928e-05 ***
## pr 1 27209 175230 289.46 4.4116 0.0463857 *
## ne 1 114114 262135 302.35 18.5023 0.0002454 ***
## ct 1 15200 163222 287.19 2.4645 0.1295338
## cum.n 1 53288 201310 293.90 8.6401 0.0071635 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1nuc.4 <- update(nuc.3, .~. - ct)
drop1(nuc.4, test="F")## Single term deletions
##
## Model:
## cost ~ date + t2 + cap + pr + ne + cum.n
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 163222 287.19
## date 1 374577 537799 323.34 57.3725 6.272e-08 ***
## t2 1 47075 210297 293.30 7.2103 0.012685 *
## cap 1 157450 320672 306.80 24.1160 4.695e-05 ***
## pr 1 42270 205492 292.56 6.4743 0.017499 *
## ne 1 127488 290709 303.66 19.5268 0.000168 ***
## cum.n 1 49676 212898 293.69 7.6087 0.010703 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(nuc.4)##
## Call:
## lm(formula = cost ~ date + t2 + cap + pr + ne + cum.n, data = nuclear)
##
## Residuals:
## Min 1Q Median 3Q Max
## -152.851 -53.929 -8.827 53.382 155.581
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.702e+03 1.294e+03 -7.495 7.55e-08 ***
## date 1.396e+02 1.843e+01 7.574 6.27e-08 ***
## t2 4.905e+00 1.827e+00 2.685 0.012685 *
## cap 4.137e-01 8.425e-02 4.911 4.70e-05 ***
## pr -8.851e+01 3.479e+01 -2.544 0.017499 *
## ne 1.502e+02 3.400e+01 4.419 0.000168 ***
## cum.n -7.919e+00 2.871e+00 -2.758 0.010703 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 80.8 on 25 degrees of freedom
## Multiple R-squared: 0.8181, Adjusted R-squared: 0.7744
## F-statistic: 18.74 on 6 and 25 DF, p-value: 3.796e-08Asemănător ca la Forward Selection, la fiecare pas se generează un nou raport cu privire la eliminarea câte unei variabile. Pentru fiecare raport, se alege p-value-ul cel mai mare și se actualizează noul model. Se repetă pașii cât timp obțin o valoare suficient de bună pentru \(p-value\) (aleasă, de obicei, cu o margine superioară de \(0.01 - 0.05\)), altfel acesta se oprește și am obținut modelul final.
Trebuie făcută observația că pentru același set de date, se obține un model diferit la Forward Selection față de Bacward Selection. Motivul principal este acela că variabilele interacționează și se influențează reciproc, de aceea coeficienții calculați își schimbă valorile în funcție de situație. De exemplu, la Forward Selection la pasul, 3 \(pt\) a fost adăugat, dar la Backward Selection, acesta fost eliminat la pasul 2. Asta înseamnă că la ultimul model, contribuția lui \(pt\) la predicție este deja adusă de alți termeni prezenți. În cadrul modelului redus, acest lucru nu se întâmpla și de aceea \(pt\) era cea mai bună opțiune la pasul respectiv.
3.2.4 Stepwise AIC Selection
Acest algoritm este o combinație a celor precedente, în sensul că la fiecare pas o variabilă poate fi adăugată sau eliminată din modelul curent. Criteriul după care are loc alegerea este calculat prin formula:
\[ AIC = -2L + 2(p + 2) \]
unde:
- \(L\) - măsură a testului de concordanță („goodness-of-fit”)
- \(p\) - numărul de parametri ai regresiei
Valoarea \(L\) si deci si \(AIC\) nu au nicio interpretare in mod direct, si se folosesc doar in compararea valorii \(AIC\) ai altui model, alegand cea mai mica valoare. După cum am menționat anterior, modelul Stepwise permite opțiunea de a șterge sau adăuga o variabilă. Acest lucru poate fi implementat utilizând funcțiile predefinite \(add1\) și \(drop1\), dar o variantă mai facilă este dată de funcția \(step\) având parametrii:
- object - obiectul cu regresia liniară deja construită
- scope - formula cu termenii care pot fi eliminați sau adăgați din regresia liniară
Exemplu: Pentru a ilustra funcționarea algoritmului, se va folosi, în continuare, setul de date mtcars ce cuprinde consumul de combustibil și 10 aspecte ale designului și performanței de automobile din anii ’70. Scopul este de a crea o regresie liniară care să prezică consumul autoturismului.
data(mtcars)
car.null <- lm(mpg ~ 1, data = mtcars) # se construiește modelul trivial (numit și modelul null)
# se definește parametrul scope ca fiind cel mai complex model, având variabile care interacționează între ele - wt, hp, cyl, disp - și cele independente -am, gear, drat, vs, qsec, carb-
car.step <- step(car.null, scope = .~. + wt * hp * factor(cyl) * disp + am
+ factor(gear) + drat + vs + qsec + carb)## Start: AIC=115.94
## mpg ~ 1
##
## Df Sum of Sq RSS AIC
## + wt 1 847.73 278.32 73.217
## + disp 1 808.89 317.16 77.397
## + factor(cyl) 2 824.78 301.26 77.752
## + hp 1 678.37 447.67 88.427
## + drat 1 522.48 603.57 97.988
## + vs 1 496.53 629.52 99.335
## + factor(gear) 2 483.24 642.80 102.003
## + am 1 405.15 720.90 103.672
## + carb 1 341.78 784.27 106.369
## + qsec 1 197.39 928.66 111.776
## <none> 1126.05 115.943
##
## Step: AIC=73.22
## mpg ~ wt
##
## Df Sum of Sq RSS AIC
## + factor(cyl) 2 95.26 183.06 63.810
## + hp 1 83.27 195.05 63.840
## + qsec 1 82.86 195.46 63.908
## + vs 1 54.23 224.09 68.283
## + carb 1 44.60 233.72 69.628
## + disp 1 31.64 246.68 71.356
## + factor(gear) 2 40.37 237.95 72.202
## <none> 278.32 73.217
## + drat 1 9.08 269.24 74.156
## + am 1 0.00 278.32 75.217
## - wt 1 847.73 1126.05 115.943
##
## Step: AIC=63.81
## mpg ~ wt + factor(cyl)
##
## Df Sum of Sq RSS AIC
## + hp 1 22.281 160.78 61.657
## + wt:factor(cyl) 2 27.170 155.89 62.669
## <none> 183.06 63.810
## + qsec 1 10.949 172.11 63.837
## + carb 1 9.244 173.81 64.152
## + vs 1 1.842 181.22 65.487
## + disp 1 0.110 182.95 65.791
## + am 1 0.090 182.97 65.794
## + drat 1 0.073 182.99 65.798
## + factor(gear) 2 6.682 176.38 66.620
## - factor(cyl) 2 95.263 278.32 73.217
## - wt 1 118.204 301.26 77.752
##
## Step: AIC=61.66
## mpg ~ wt + factor(cyl) + hp
##
## Df Sum of Sq RSS AIC
## + wt:hp 1 34.623 126.16 55.897
## + hp:factor(cyl) 2 28.550 132.23 59.401
## + wt:factor(cyl) 2 25.385 135.39 60.158
## + am 1 9.752 151.03 61.655
## <none> 160.78 61.657
## + drat 1 2.438 158.34 63.168
## + carb 1 1.262 159.51 63.405
## + vs 1 0.655 160.12 63.527
## + disp 1 0.651 160.13 63.527
## + qsec 1 0.229 160.55 63.611
## - hp 1 22.281 183.06 63.810
## - factor(cyl) 2 34.270 195.05 63.840
## + factor(gear) 2 7.366 153.41 64.156
## - wt 1 116.390 277.17 77.084
##
## Step: AIC=55.9
## mpg ~ wt + factor(cyl) + hp + wt:hp
##
## Df Sum of Sq RSS AIC
## - factor(cyl) 2 3.606 129.76 52.799
## <none> 126.16 55.897
## + qsec 1 5.278 120.88 56.529
## + vs 1 2.345 123.81 57.297
## + hp:factor(cyl) 2 8.645 117.51 57.625
## + drat 1 0.229 125.93 57.839
## + am 1 0.082 126.07 57.876
## + carb 1 0.038 126.12 57.887
## + disp 1 0.032 126.12 57.889
## + factor(gear) 2 7.376 118.78 57.969
## + wt:factor(cyl) 2 1.070 125.08 59.624
## - wt:hp 1 34.623 160.78 61.657
##
## Step: AIC=52.8
## mpg ~ wt + hp + wt:hp
##
## Df Sum of Sq RSS AIC
## + qsec 1 8.720 121.04 52.573
## <none> 129.76 52.799
## + vs 1 4.324 125.44 53.714
## + disp 1 0.591 129.17 54.653
## + carb 1 0.176 129.59 54.755
## + am 1 0.042 129.72 54.788
## + drat 1 0.000 129.76 54.799
## + factor(gear) 2 7.219 122.54 54.967
## + factor(cyl) 2 3.606 126.16 55.897
## - wt:hp 1 65.286 195.05 63.840
##
## Step: AIC=52.57
## mpg ~ wt + hp + qsec + wt:hp
##
## Df Sum of Sq RSS AIC
## <none> 121.04 52.573
## - qsec 1 8.720 129.76 52.799
## + factor(gear) 2 9.482 111.56 53.962
## + am 1 1.939 119.10 54.056
## + carb 1 0.080 120.96 54.551
## + drat 1 0.012 121.03 54.570
## + vs 1 0.010 121.03 54.570
## + disp 1 0.008 121.03 54.571
## + factor(cyl) 2 0.164 120.88 56.529
## - wt:hp 1 65.018 186.06 64.331summary(car.step)##
## Call:
## lm(formula = mpg ~ wt + hp + qsec + wt:hp, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.8243 -1.3980 0.0303 1.1582 4.3650
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 40.310410 7.677887 5.250 1.56e-05 ***
## wt -8.681516 1.292525 -6.717 3.28e-07 ***
## hp -0.106181 0.026263 -4.043 0.000395 ***
## qsec 0.503163 0.360768 1.395 0.174476
## wt:hp 0.027791 0.007298 3.808 0.000733 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.117 on 27 degrees of freedom
## Multiple R-squared: 0.8925, Adjusted R-squared: 0.8766
## F-statistic: 56.05 on 4 and 27 DF, p-value: 1.094e-12Funcția \(step\) returnează un model antrenat și afișează un raport detaliat în legătură cu toți pașii făcuți. La fiecare pas, algoritmul alege cea mai bună variantă dpdv al coeficientului AIC și adaugă/elimină variabila corespunzătoare. Acesta se sfârșește atunci când cea mai bună alegere este să nu se mai modifice termenii regresiei.
Forula regresiei liniare la care a ajuns algoritmul este: \(mpg \sim wt + hp + qsec + wt:hp\)
3.3 Interpretarea datelor și observațiilor
Este important să se investigheze orice observație care pare neobișnuită sau extremă în comparație cu restul setului de date. În general se recomandă o analiză exploratorie a datelor, care implică statistici, grafice, tabele deoarece se identifică mai ușor astfel de valori extreme care pot afecta în mod negativ cercetările ulterioare. Astfel, se definesc termeni importanți precum:
- Valoare aberantă (outlier) - observație care arată neobișnuit în contextul unui set de date. De obicei, în cazul regresiilor liniare, o valoare aberantă are o valoare reziduală mare, dar cu condiția ca punctul să nu se conformeze tendinței modelului
- Punct pârghie (leverage) - o observație a cărei valoare x este neobișnuită, dar a cărei valoare prezisă y se conformează cu trendul și urmează linia de regresie prezisă.
- Valoare influentă (influence) - o observație cu un efect de levier mare care afectează tendințele estimate
Exemple care ilustrează termenii definiți:
x <- c(1.1, 1.3, 2.3, 1.6, 1.2, 0.1, 1.8, 1.9, 0.2, 0.75)
y <- c(6.7, 7.9, 9.8, 9.3, 8.2, 2.9, 6.6, 11.1, 4.7, 3)
p1x <- 1.2
p1y <- 14
p2x <- 5
p2y <- 19
p3x <- 5
p3y <- 5
mod.0 <- lm(y~x)
mod.1 <- lm(c(y, p1y) ~ c(x, p1x))
mod.2 <- lm(c(y, p2y) ~ c(x, p2x))
mod.3 <- lm(c(y, p3y) ~ c(x, p3x))
plot(x,y,xlim=c(0,5),ylim=c(0,20))
points(p1x,p1y,pch=15,cex=1.5)
abline(mod.0)
abline(mod.1,lty=2)
text(2.4,1,labels="Valoare aberanta(outlier), efect de levier(leverage) mic, influenta mica",cex=1.1)
plot(x,y,xlim=c(0,5),ylim=c(0,20))
points(p2x,p2y,pch=15,cex=1.5)
abline(mod.0)
abline(mod.2,lty=2)
text(2.4,1,labels="Nu este valoare aberanta(not outlier), efect de levier(leverage) mare, influenta mica",cex=0.9)
plot(x,y,xlim=c(0,5),ylim=c(0,20))
points(p3x,p3y,pch=15,cex=1.5)
abline(mod.0)
abline(mod.3,lty=2)
text(2.4,1,labels="Valoare aberanta(outlier), efect de levier(leverage) mare, influenta mare",cex=1.1)
Punctul îngroșat reprezintă observația adăugată care este clasificată (valoare aberantă, punct pârghie, punct influent). Linia punctată reprezintă noua dreapta de regresie obținută în urma adăugării ultimei observații.
- În prima figură, observația prezintă o valoare aberantă mare deoarece predicția nu se conformează cu dreapta de regresie. Cu toate acestea, are un efect de levier mic (low leverage) datorită valorii sale y care este asemăntoare cu restul observațiilor
- În figura 2 observația nu este o valoare aberantă deoarece se conformează cu dreapte de regresie liniară și este un punct de pârghie deoarece are o valoare de predicție extremă (mult mai mare față de restul datelor). Se observă ca nu influențează puternic dreapta de regresie liniară
- În figura 3, se observă că observația reprezintă o valoare aberantă din moment ce are o valoare reziduală mare și nu se conformează cu dreapta de regresie liniară. În plus, este un punct pârghie deoarece are o valoare de predicție extremă. Observația influențează modelul de regresie liniară, schimbarea direcției si pantei fiind evidentă
3.3.1 Puncte pârghie (leverage)
Această metrică se calculează folosind matricea de structură definită în capitolele trecute. Concret, pentru cele \(n\) observații, efectul de levier pentru observația \(i\) este notat \(h_{ii}\). Valorile pentru cele \(n\) observațtii sunt cele de pe diagonala principală a matricei:
\[ H = X(X^TX)^{-1}X^T \]
leverage <- hatvalues(mod.0)
plot(leverage)
3.3.2 Distanța Cook
O altă metrică importantă pentru determinarea gradului de influență a unei observații o reprezintă distanța Cook care estimează amploarea efectului de eliminare a observației \(i\) din modelul adaptat.
Distanța Cook pentru observația \(D_{i}\) (masură a influenței celei de-a \(i\)-a observații asupra tutoror valorilor prognozate) se calculează:
\[D_{i} = \sum_{j=1}^{n} \frac{(\hat{y_{j}} - \hat{y_{j}}^{(-i)})^2}{(p+1)\hat{\sigma}^{2}}, \space \space \space \space i,j=1, \dots ,n\]
- \(\hat{\sigma}^{2}\) - deviația standard a erorii
- \(\hat{y_{j}}\) - valoarea estimată pentru \(j\)-a observație
- \(\hat{y_{j}}^{(-i)}\) - valoarea estimată din regresia calculată după omiterea cele de-a \(i\)-a observații
- \(p\) - numărul de parametrii ai regresiei
- \(n\) - numărul de observații
Pentru a considera că a \(i\)-a observație este suficient de influentă pentru modelul de regresie, pentru se ia o valoare de prag egală, în general, cu \(4/n\). Astfel, dacă \(D_{i} > 4/n\), atunci observația \(i\) se poate considera suficient de influentă.
Exemplu: Pentru cele 3 modele de regresie definite anterior \(mod1\), \(mod2\), \(mod3\) calculăm distanța Cook folosind funcția predefintă din R:
n <- length(mod.1)
plot(mod.1,which=4)
abline(h=c(1,4/n),lty=2)
plot(mod.2,which=4)
abline(h=c(1,4/n),lty=2)
plot(mod.3,which=4)
abline(h=c(1,4/n),lty=2)
Se observă că în figura a 3-a, distanța Cook a ultimei observații depășește pragul de \(4/n\), deci se poate considera că observația respectivă influențează modelul de regresie.
3.4 Combinarea metricilor prezentate
Graficele prezentate până acum prezintă statistici cu privire la valorea reziduală standardizată, efectul de levier (leverage) și distanța Cook pentru cea de-a \(i\)-a observație. Pentru o mai bună înțelegere a observațiilor și efectele pe care acestea le au asupra regresiilor liniare, aceste metrici pot fi combinate într-o singură reprezentare grafică.
Se folosește în continuare setul de date definit anterior \(mod.1\), \(mod.2\), \(mod.3\). Pentru fiecare grafic, se afișează efectul de pârghie (leverage) pe axa orizontală și reziduurile standardizate pe axa verticală pentru fiecare observație. Fiind o funcție a reziduurilor si a efectului de pârghie, distanțele Cook pot fi reprezentate ca și contururi. Aceste contururi delimitează zonele spațiale ale porțiunilor care corespund unei influențe mari (colțurile extreme din dreapta).
plot(mod.1, which = 5, add.smooth = FALSE, cook.levels = c(4/11, 0.5, 1))
plot(mod.2, which = 5, add.smooth = FALSE, cook.levels = c(4/11, 0.5, 1))
plot(mod.3, which = 5, add.smooth = FALSE, cook.levels = c(4/11, 0.5, 1))
Cu cât un punct se apropie de linia orizontală, cu atât valoarea reziduală este mai mică. Un punct care se află mult mai la stânga decât la dreapta are un efect de levier (leverage) mai mic. Dacă un punct se află suficient de departe de linia orizontală, având în vedere poziția sa de pârghie (axa x), acesta va depăși contururile care marchează anumite valori ale lui \(D_{i}\), indicând o influență mare.
Graficul pentru \(mod.1\) arată că punctul adăugat are un reziduu mare, dar nu incalcă conturul \(4/11\) deoarece se află într-o poziție de levier(leverage) scăzut. Graficul pentru \(mod.2\) arată că punctul adăugat se află într-o poziție de levier(leverage) ridicat, dar nu este influent deoarece reziduul său este mic. În cele din urmă, grafiul pentru \(mod.3\) identifică în mod clar punctul adăugat ca fiind extrem de influent - cu o valoare reziduală mare și efect de levier (leverage) ridicat.
Un alt tip de grafic care prezintă aceleași informații, doar că pe axa verticală este afișată \(distanța \space Cook\), iar pe cea orizontală o transformarea a efectului levier (leverage) dată de formula \(h_{ii}/(1-h_{ii})\). Această transformare amplifică punctele de pârghie mai mari în ceea ce privește poziția lor orizontală. Liniile punctate reprezintă reziduurile standardizate (ca o funcție dintre efectul de levier(leverage) și reziduuri) și etichetate cu magnitudinea lor.
plot(mod.1, which = 6, add.smooth = FALSE)
abline(h = c(4 / 11), lty = 1, col = 'red')
plot(mod.2, which = 6, add.smooth = FALSE)
abline(h = c(4 / 11), lty = 1, col = 'red')
plot(mod.3, which = 6, add.smooth = FALSE)
abline(h = c(4 / 11), lty = 1, col = 'red')
În figuri a fost adăugată și dreapta pentru distanța Cook aferentă valorii de prag \(D_{i} = 4/11\). Astfel, pentru fiecare punct care se află deasupra liniei respective, putem considera că acesta are un grad de influență ridicat pentru modelul de regresie liniară aferent.
Exercițiul 20.1
a.
#library("boot")
nuc.null <- lm(cost ~ 1,data = nuclear)
nuc.step <- step(nuc.null, scope= .~. + date +t1 + t2 + cap + pr + ne + ct + bw + cum.n + pt)## Start: AIC=329.72
## cost ~ 1
##
## Df Sum of Sq RSS AIC
## + date 1 334335 562837 316.80
## + pt 1 305334 591839 318.41
## + cap 1 199673 697499 323.66
## + t1 1 186984 710189 324.24
## + ne 1 128641 768531 326.77
## + cum.n 1 67938 829234 329.20
## <none> 897172 329.72
## + ct 1 43042 854130 330.15
## + bw 1 16205 880967 331.14
## + pr 1 9037 888136 331.40
## + t2 1 27 897145 331.72
##
## Step: AIC=316.8
## cost ~ date
##
## Df Sum of Sq RSS AIC
## + cap 1 189732 373105 305.64
## + ne 1 92256 470581 313.07
## + pt 1 90587 472250 313.18
## + t2 1 68161 494676 314.67
## + ct 1 54794 508043 315.52
## <none> 562837 316.80
## + t1 1 15322 547515 317.92
## + cum.n 1 4658 558179 318.53
## + pr 1 4027 558810 318.57
## + bw 1 1240 561597 318.73
## - date 1 334335 897172 329.72
##
## Step: AIC=305.64
## cost ~ date + cap
##
## Df Sum of Sq RSS AIC
## + pt 1 95917 277189 298.13
## + ne 1 94537 278569 298.29
## + ct 1 48957 324148 303.14
## + cum.n 1 28062 345044 305.14
## <none> 373105 305.64
## + pr 1 18453 354652 306.02
## + t2 1 13356 359750 306.48
## + bw 1 7505 365601 306.99
## + t1 1 957 372148 307.56
## - cap 1 189732 562837 316.80
## - date 1 324394 697499 323.66
##
## Step: AIC=298.13
## cost ~ date + cap + pt
##
## Df Sum of Sq RSS AIC
## + ne 1 54572 222617 293.12
## + ct 1 18119 259069 297.97
## <none> 277189 298.13
## + t2 1 10429 266760 298.91
## + pr 1 6017 271171 299.43
## + bw 1 654 276535 300.06
## + cum.n 1 463 276726 300.08
## + t1 1 12 277177 300.13
## - pt 1 95917 373105 305.64
## - date 1 111397 388586 306.94
## - cap 1 195062 472250 313.19
##
## Step: AIC=293.12
## cost ~ date + cap + pt + ne
##
## Df Sum of Sq RSS AIC
## + t2 1 19230 203387 292.23
## + ct 1 15765 206852 292.77
## + cum.n 1 13820 208797 293.07
## <none> 222617 293.12
## + pr 1 5231 217386 294.36
## + bw 1 448 222169 295.06
## + t1 1 107 222510 295.10
## - ne 1 54572 277189 298.13
## - pt 1 55952 278569 298.29
## - date 1 119333 341950 304.85
## - cap 1 195756 418373 311.31
##
## Step: AIC=292.23
## cost ~ date + cap + pt + ne + t2
##
## Df Sum of Sq RSS AIC
## + pr 1 23244 180143 290.35
## + cum.n 1 12951 190436 292.12
## <none> 203387 292.23
## + ct 1 10212 193175 292.58
## - t2 1 19230 222617 293.12
## + bw 1 860 202527 294.09
## + t1 1 368 203019 294.17
## - pt 1 50299 253686 297.30
## - ne 1 63373 266760 298.91
## - cap 1 132791 336178 306.31
## - date 1 138352 341739 306.83
##
## Step: AIC=290.34
## cost ~ date + cap + pt + ne + t2 + pr
##
## Df Sum of Sq RSS AIC
## + cum.n 1 20987 159156 288.38
## <none> 180143 290.35
## + t1 1 6495 173648 291.17
## + bw 1 5898 174245 291.28
## + ct 1 4828 175316 291.48
## - pr 1 23244 203387 292.23
## - pt 1 32754 212898 293.69
## - t2 1 37243 217386 294.36
## - ne 1 67232 247375 298.49
## - cap 1 131615 311758 305.90
## - date 1 158588 338731 308.55
##
## Step: AIC=288.38
## cost ~ date + cap + pt + ne + t2 + pr + cum.n
##
## Df Sum of Sq RSS AIC
## - pt 1 4066 163222 287.19
## + ct 1 11794 147362 287.92
## <none> 159156 288.38
## + t1 1 3103 156053 289.75
## + bw 1 2806 156350 289.81
## - cum.n 1 20987 180143 290.35
## - pr 1 31280 190436 292.12
## - t2 1 40640 199796 293.66
## - ne 1 86979 246135 300.33
## - date 1 150115 309271 307.64
## - cap 1 150213 309369 307.65
##
## Step: AIC=287.19
## cost ~ date + cap + ne + t2 + pr + cum.n
##
## Df Sum of Sq RSS AIC
## + ct 1 15200 148021 286.06
## <none> 163222 287.19
## + bw 1 4884 158337 288.22
## + pt 1 4066 159156 288.38
## + t1 1 2452 160770 288.70
## - pr 1 42270 205492 292.56
## - t2 1 47075 210297 293.30
## - cum.n 1 49676 212898 293.69
## - ne 1 127488 290709 303.66
## - cap 1 157450 320672 306.80
## - date 1 374577 537799 323.34
##
## Step: AIC=286.06
## cost ~ date + cap + ne + t2 + pr + cum.n + ct
##
## Df Sum of Sq RSS AIC
## <none> 148021 286.06
## - ct 1 15200 163222 287.19
## + t1 1 2853 145168 287.44
## + bw 1 1661 146360 287.70
## + pt 1 660 147362 287.92
## - pr 1 27209 175230 289.46
## - t2 1 28539 176560 289.70
## - cum.n 1 53288 201310 293.90
## - ne 1 114114 262135 302.35
## - cap 1 162846 310867 307.80
## - date 1 369019 517040 324.08summary(nuc.step)##
## Call:
## lm(formula = cost ~ date + cap + ne + t2 + pr + cum.n + ct, data = nuclear)
##
## Residuals:
## Min 1Q Median 3Q Max
## -137.833 -46.897 -3.834 38.100 178.174
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.606e+03 1.259e+03 -7.627 7.27e-08 ***
## date 1.386e+02 1.792e+01 7.735 5.70e-08 ***
## cap 4.215e-01 8.203e-02 5.138 2.93e-05 ***
## ne 1.434e+02 3.333e+01 4.301 0.000245 ***
## t2 4.011e+00 1.865e+00 2.151 0.041748 *
## pr -7.372e+01 3.510e+01 -2.100 0.046386 *
## cum.n -8.222e+00 2.797e+00 -2.939 0.007164 **
## ct 4.757e+01 3.030e+01 1.570 0.129534
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 78.53 on 24 degrees of freedom
## Multiple R-squared: 0.835, Adjusted R-squared: 0.7869
## F-statistic: 17.35 on 7 and 24 DF, p-value: 5.663e-08b. Modelele difera.
Forward Selection: lm(formula = cost ~ date + cap + pt + ne, data = nuclear)
Backward Selection: lm(formula = cost ~ date + t2 + cap + pr + ne + cum.n, data = nuclear)
AIC Selection: lm(formula = cost ~ date + cap + ne + t2 + pr + cum.n + ct, data = nuclear)
Se observă ca diferă prin mai multe variabile față de Forward Selection, iar față de Backward Selection doar prin una (\(ct\))
c.
gal <- data.frame(d = c(573, 534, 495, 451, 395, 337, 253), h = c(1000, 800., 600, 450, 300, 200, 100))
gal.0 <- lm(d ~ 1, data = gal)
gal.1 <- lm(d ~ h, data = gal)
gal.2 <- lm(d ~ h + I(h^2), data = gal)
gal.3 <- lm(d ~ h + I(h^2) + I(h^3), data = gal)
gal.4 <- lm(d ~ h + I(h^2) + I(h^3) + I(h^4), data = gal)d.
anova(gal.0, gal.1, gal.2, gal.3, gal.4)## Analysis of Variance Table
##
## Model 1: d ~ 1
## Model 2: d ~ h
## Model 3: d ~ h + I(h^2)
## Model 4: d ~ h + I(h^2) + I(h^3)
## Model 5: d ~ h + I(h^2) + I(h^3) + I(h^4)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 6 77022
## 2 5 5671 1 71351 11208.1079 8.921e-05 ***
## 3 4 744 1 4927 773.9758 0.001290 **
## 4 3 48 1 696 109.3033 0.009025 **
## 5 2 13 1 36 5.5799 0.142011
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1e)
#library("faraway")
data("diabetes")
diab.filtered <- na.omit(diabetes[,c("chol", "age", "gender", "height", "weight", "frame", "waist", "hip", "location")])f.
dia.null <- lm(chol ~ 1, data = diab.filtered)
dia.full <- lm(chol ~ age * gender * weight * frame + waist * height * hip + location, data = diab.filtered)g.
dia.aic <- step(dia.null, scope = .~. + age * gender * weight * frame + waist * height * hip + location)## Start: AIC=2892.63
## chol ~ 1
##
## Df Sum of Sq RSS AIC
## + age 1 36863 689181 2874.7
## + frame 2 20262 705781 2885.8
## + waist 1 13442 712601 2887.5
## + hip 1 5567 720476 2891.7
## + weight 1 4954 721090 2892.0
## <none> 726044 2892.6
## + location 1 2511 723532 2893.3
## + height 1 1891 724153 2893.6
## + gender 1 494 725550 2894.4
##
## Step: AIC=2874.67
## chol ~ age
##
## Df Sum of Sq RSS AIC
## + frame 2 17740 671441 2868.7
## + waist 1 7079 682102 2872.7
## + weight 1 6147 683034 2873.2
## + hip 1 4774 684408 2874.0
## <none> 689181 2874.7
## + location 1 2676 686505 2875.2
## + gender 1 1089 688092 2876.1
## + height 1 660 688521 2876.3
## - age 1 36863 726044 2892.6
##
## Step: AIC=2868.68
## chol ~ age + frame
##
## Df Sum of Sq RSS AIC
## + age:frame 2 9896 661545 2867.0
## + waist 1 5820 665621 2867.3
## + weight 1 5007 666434 2867.8
## <none> 671441 2868.7
## + hip 1 2780 668661 2869.1
## + location 1 1911 669530 2869.6
## + gender 1 736 670705 2870.3
## + height 1 282 671159 2870.5
## - frame 2 17740 689181 2874.7
## - age 1 34341 705781 2885.8
##
## Step: AIC=2867
## chol ~ age + frame + age:frame
##
## Df Sum of Sq RSS AIC
## + waist 1 4468.8 657076 2866.4
## <none> 661545 2867.0
## + weight 1 3207.1 658338 2867.1
## + hip 1 2022.8 659522 2867.8
## + location 1 1660.4 659884 2868.0
## + gender 1 623.4 660921 2868.6
## - age:frame 2 9896.1 671441 2868.7
## + height 1 355.2 661189 2868.8
##
## Step: AIC=2866.4
## chol ~ age + frame + waist + age:frame
##
## Df Sum of Sq RSS AIC
## <none> 657076 2866.4
## - waist 1 4468.8 661545 2867.0
## - age:frame 2 8545.1 665621 2867.3
## + location 1 1585.1 655491 2867.5
## + height 1 442.5 656633 2868.1
## + hip 1 262.7 656813 2868.2
## + gender 1 225.4 656850 2868.3
## + weight 1 3.6 657072 2868.4summary(dia.aic)##
## Call:
## lm(formula = chol ~ age + frame + waist + age:frame, data = diab.filtered)
##
## Residuals:
## Min 1Q Median 3Q Max
## -126.228 -26.008 -5.256 21.263 221.888
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 131.3260 18.3762 7.147 4.66e-12 ***
## age 0.9939 0.2662 3.734 0.000218 ***
## framemedium 28.8524 15.4951 1.862 0.063377 .
## framelarge 39.1384 19.4922 2.008 0.045369 *
## waist 0.6875 0.4299 1.599 0.110633
## age:framemedium -0.3788 0.3336 -1.136 0.256864
## age:framelarge -0.8286 0.3754 -2.208 0.027880 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 41.8 on 376 degrees of freedom
## Multiple R-squared: 0.09499, Adjusted R-squared: 0.08055
## F-statistic: 6.578 on 6 and 376 DF, p-value: 1.271e-06h.
add1(dia.null, scope = .~. + age * gender * weight * frame + waist * height * hip + location, test = "F")## Single term additions
##
## Model:
## chol ~ 1
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 726044 2892.6
## age 1 36863 689181 2874.7 20.3787 8.48e-06 ***
## gender 1 494 725550 2894.4 0.2593 0.610892
## weight 1 4954 721090 2892.0 2.6175 0.106518
## frame 2 20262 705781 2885.8 5.4547 0.004618 **
## waist 1 13442 712601 2887.5 7.1871 0.007662 **
## height 1 1891 724153 2893.6 0.9949 0.319187
## hip 1 5567 720476 2891.7 2.9441 0.087004 .
## location 1 2511 723532 2893.3 1.3224 0.250874
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1dia.1 <- update(dia.null, .~. + age)
add1(dia.1, scope = .~. + age * gender * weight * frame + waist * height * hip + location, test = "F")## Single term additions
##
## Model:
## chol ~ age
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 689181 2874.7
## gender 1 1089.1 688092 2876.1 0.6015 0.438497
## weight 1 6147.4 683034 2873.2 3.4200 0.065186 .
## frame 2 17740.3 671441 2868.7 5.0068 0.007142 **
## waist 1 7079.1 682102 2872.7 3.9438 0.047763 *
## height 1 659.6 688521 2876.3 0.3640 0.546629
## hip 1 4773.6 684408 2874.0 2.6504 0.104353
## location 1 2676.4 686505 2875.2 1.4815 0.224302
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1dia.2 <- update(dia.1, .~. + frame)
add1(dia.2, scope = .~. + age * gender * weight * frame + waist * height * hip + location, test = "F")## Single term additions
##
## Model:
## chol ~ age + frame
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 671441 2868.7
## gender 1 735.5 670705 2870.3 0.4145 0.52007
## weight 1 5007.3 666434 2867.8 2.8401 0.09276 .
## waist 1 5819.8 665621 2867.3 3.3050 0.06986 .
## height 1 281.6 671159 2870.5 0.1586 0.69068
## hip 1 2779.6 668661 2869.1 1.5714 0.21078
## location 1 1910.8 669530 2869.6 1.0788 0.29964
## age:frame 2 9896.1 661545 2867.0 2.8198 0.06088 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Având \(alpha = 0.05\), algoritmul se oprește.
Modelul obținut la subpunctul g): chol ~ age + frame + waist + age:frame
Modelul obținut la subpunctul h): chol ~ age + frame
Se observă că Stepwise AIC a adăugat în plus factorul \(waist\) și interacțiunea dintre \(age\) și \(frame\). Forward Selection a considerat că nu se obține o îmbunătățire suficient de bună astfel încât acestea să fie luate în considerare
i.
dia.aic2 <- step(dia.full)## Start: AIC=2876.45
## chol ~ age * gender * weight * frame + waist * height * hip +
## location
##
## Df Sum of Sq RSS AIC
## - age:gender:weight:frame 2 1570.6 593562 2873.5
## - location 1 1704.8 593696 2875.6
## <none> 591991 2876.4
## - waist:height:hip 1 4221.1 596212 2877.2
##
## Step: AIC=2873.47
## chol ~ age + gender + weight + frame + waist + height + hip +
## location + age:gender + age:weight + gender:weight + age:frame +
## gender:frame + weight:frame + waist:height + waist:hip +
## height:hip + age:gender:weight + age:gender:frame + age:weight:frame +
## gender:weight:frame + waist:height:hip
##
## Df Sum of Sq RSS AIC
## - age:gender:frame 2 1473.5 595035 2870.4
## - age:weight:frame 2 4008.6 597570 2872.0
## - gender:weight:frame 2 4114.0 597676 2872.1
## - location 1 1448.8 595011 2872.4
## - age:gender:weight 1 1883.2 595445 2872.7
## <none> 593562 2873.5
## - waist:height:hip 1 4131.1 597693 2874.1
##
## Step: AIC=2870.41
## chol ~ age + gender + weight + frame + waist + height + hip +
## location + age:gender + age:weight + gender:weight + age:frame +
## gender:frame + weight:frame + waist:height + waist:hip +
## height:hip + age:gender:weight + age:weight:frame + gender:weight:frame +
## waist:height:hip
##
## Df Sum of Sq RSS AIC
## - age:weight:frame 2 4209.1 599244 2869.1
## - location 1 1482.2 596517 2869.4
## - gender:weight:frame 2 5363.1 600398 2869.8
## - age:gender:weight 1 2751.8 597787 2870.2
## <none> 595035 2870.4
## - waist:height:hip 1 3989.8 599025 2871.0
##
## Step: AIC=2869.11
## chol ~ age + gender + weight + frame + waist + height + hip +
## location + age:gender + age:weight + gender:weight + age:frame +
## gender:frame + weight:frame + waist:height + waist:hip +
## height:hip + age:gender:weight + gender:weight:frame + waist:height:hip
##
## Df Sum of Sq RSS AIC
## - age:frame 2 386.7 599631 2865.4
## - location 1 1258.2 600503 2867.9
## - gender:weight:frame 2 5609.4 604854 2868.7
## - age:gender:weight 1 2882.2 602127 2868.9
## <none> 599244 2869.1
## - waist:height:hip 1 3980.4 603225 2869.7
##
## Step: AIC=2865.36
## chol ~ age + gender + weight + frame + waist + height + hip +
## location + age:gender + age:weight + gender:weight + gender:frame +
## weight:frame + waist:height + waist:hip + height:hip + age:gender:weight +
## gender:weight:frame + waist:height:hip
##
## Df Sum of Sq RSS AIC
## - location 1 1284.0 600915 2864.2
## - gender:weight:frame 2 5824.2 605455 2865.1
## - age:gender:weight 1 3030.1 602661 2865.3
## <none> 599631 2865.4
## - waist:height:hip 1 4146.7 603778 2866.0
##
## Step: AIC=2864.18
## chol ~ age + gender + weight + frame + waist + height + hip +
## age:gender + age:weight + gender:weight + gender:frame +
## weight:frame + waist:height + waist:hip + height:hip + age:gender:weight +
## gender:weight:frame + waist:height:hip
##
## Df Sum of Sq RSS AIC
## - gender:weight:frame 2 5701.2 606616 2863.8
## - age:gender:weight 1 2977.2 603892 2864.1
## <none> 600915 2864.2
## - waist:height:hip 1 3945.4 604860 2864.7
##
## Step: AIC=2863.8
## chol ~ age + gender + weight + frame + waist + height + hip +
## age:gender + age:weight + gender:weight + gender:frame +
## weight:frame + waist:height + waist:hip + height:hip + age:gender:weight +
## waist:height:hip
##
## Df Sum of Sq RSS AIC
## - weight:frame 2 3371.5 609988 2861.9
## - gender:frame 2 4023.3 610640 2862.3
## - waist:height:hip 1 1577.2 608193 2862.8
## <none> 606616 2863.8
## - age:gender:weight 1 3460.2 610077 2864.0
##
## Step: AIC=2861.92
## chol ~ age + gender + weight + frame + waist + height + hip +
## age:gender + age:weight + gender:weight + gender:frame +
## waist:height + waist:hip + height:hip + age:gender:weight +
## waist:height:hip
##
## Df Sum of Sq RSS AIC
## - gender:frame 2 3209.7 613197 2859.9
## - waist:height:hip 1 1617.6 611605 2860.9
## <none> 609988 2861.9
## - age:gender:weight 1 4315.1 614303 2862.6
##
## Step: AIC=2859.93
## chol ~ age + gender + weight + frame + waist + height + hip +
## age:gender + age:weight + gender:weight + waist:height +
## waist:hip + height:hip + age:gender:weight + waist:height:hip
##
## Df Sum of Sq RSS AIC
## - waist:height:hip 1 1841.4 615039 2859.1
## <none> 613197 2859.9
## - age:gender:weight 1 4169.7 617367 2860.5
## - frame 2 7774.5 620972 2860.8
##
## Step: AIC=2859.08
## chol ~ age + gender + weight + frame + waist + height + hip +
## age:gender + age:weight + gender:weight + waist:height +
## waist:hip + height:hip + age:gender:weight
##
## Df Sum of Sq RSS AIC
## - waist:height 1 119.3 615158 2857.2
## - height:hip 1 196.6 615235 2857.2
## <none> 615039 2859.1
## - age:gender:weight 1 4835.4 619874 2860.1
## - frame 2 8260.9 623300 2860.2
## - waist:hip 1 8874.0 623913 2862.6
##
## Step: AIC=2857.15
## chol ~ age + gender + weight + frame + waist + height + hip +
## age:gender + age:weight + gender:weight + waist:hip + height:hip +
## age:gender:weight
##
## Df Sum of Sq RSS AIC
## - height:hip 1 1248.5 616407 2855.9
## <none> 615158 2857.2
## - age:gender:weight 1 4828.3 619986 2858.2
## - frame 2 8801.8 623960 2858.6
## - waist:hip 1 8950.0 624108 2860.7
##
## Step: AIC=2855.93
## chol ~ age + gender + weight + frame + waist + height + hip +
## age:gender + age:weight + gender:weight + waist:hip + age:gender:weight
##
## Df Sum of Sq RSS AIC
## - height 1 727.1 617134 2854.4
## <none> 616407 2855.9
## - age:gender:weight 1 4713.5 621120 2856.8
## - frame 2 9157.6 625564 2857.6
## - waist:hip 1 8105.1 624512 2858.9
##
## Step: AIC=2854.38
## chol ~ age + gender + weight + frame + waist + hip + age:gender +
## age:weight + gender:weight + waist:hip + age:gender:weight
##
## Df Sum of Sq RSS AIC
## <none> 617134 2854.4
## - age:gender:weight 1 4969.6 622103 2855.4
## - frame 2 9423.4 626557 2856.2
## - waist:hip 1 7822.5 624956 2857.2summary(dia.aic2)##
## Call:
## lm(formula = chol ~ age + gender + weight + frame + waist + hip +
## age:gender + age:weight + gender:weight + waist:hip + age:gender:weight,
## data = diab.filtered)
##
## Residuals:
## Min 1Q Median 3Q Max
## -129.934 -24.257 -3.982 22.115 212.604
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -90.024818 111.236018 -0.809 0.4189
## age 1.851810 1.056141 1.753 0.0804 .
## genderfemale 32.626630 60.751227 0.537 0.5916
## weight 0.549625 0.290229 1.894 0.0590 .
## framemedium 10.484266 5.542926 1.891 0.0593 .
## framelarge 0.309343 7.117136 0.043 0.9654
## waist 6.213040 2.521085 2.464 0.0142 *
## hip 3.788484 2.441865 1.551 0.1216
## age:genderfemale -0.980576 1.292190 -0.759 0.4484
## age:weight -0.011434 0.005813 -1.967 0.0499 *
## genderfemale:weight -0.474328 0.335521 -1.414 0.1583
## waist:hip -0.119344 0.055108 -2.166 0.0310 *
## age:genderfemale:weight 0.012477 0.007228 1.726 0.0852 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 40.84 on 370 degrees of freedom
## Multiple R-squared: 0.15, Adjusted R-squared: 0.1224
## F-statistic: 5.441 on 12 and 370 DF, p-value: 1.525e-08Modelul obținut este: chol ~ age + gender + weight + frame + waist + hip + age:gender + age:weight + gender:weight + waist:hip + age:gender:weight Se observă diferența, acesta fiind un model cu mult mai complex față de cel de la punctul precedent. Având ca punct de plecare un model mult mai complex față de cel trivial, algoritmul Stepwise AIC a găsit combinații și calculat valori AIC mai mici având punctul de plecare \(dia.full\) față de punctul de plecare \(dia.null\)
j.
car.1 <- lm(I(1/mpg) ~ 1, data = mtcars)
car.step <- step(car.1, scope = .~. + wt * hp * factor(cyl) * disp + am + factor(gear) + drat + vs + qsec + carb)## Start: AIC=-261.99
## I(1/mpg) ~ 1
##
## Df Sum of Sq RSS AIC
## + wt 1 0.0066224 0.0017402 -310.22
## + disp 1 0.0064733 0.0018892 -307.60
## + factor(cyl) 2 0.0055711 0.0027914 -293.10
## + hp 1 0.0048677 0.0034948 -287.91
## + vs 1 0.0034243 0.0049382 -276.85
## + drat 1 0.0034045 0.0049580 -276.72
## + factor(gear) 2 0.0034781 0.0048844 -275.20
## + am 1 0.0024375 0.0059251 -271.02
## + carb 1 0.0023167 0.0060458 -270.37
## + qsec 1 0.0012450 0.0071175 -265.15
## <none> 0.0083625 -261.99
##
## Step: AIC=-310.22
## I(1/mpg) ~ wt
##
## Df Sum of Sq RSS AIC
## + hp 1 0.0004614 0.0012787 -318.08
## + qsec 1 0.0004578 0.0012824 -317.99
## + disp 1 0.0003175 0.0014226 -314.67
## + vs 1 0.0002579 0.0014823 -313.36
## + factor(cyl) 2 0.0003238 0.0014163 -312.81
## + carb 1 0.0002176 0.0015226 -312.50
## + factor(gear) 2 0.0002346 0.0015055 -310.86
## <none> 0.0017402 -310.22
## + am 1 0.0000937 0.0016465 -310.00
## + drat 1 0.0000003 0.0017399 -308.23
## - wt 1 0.0066224 0.0083625 -261.99
##
## Step: AIC=-318.08
## I(1/mpg) ~ wt + hp
##
## Df Sum of Sq RSS AIC
## <none> 0.0012787 -318.08
## + qsec 1 0.00004906 0.0012297 -317.34
## + disp 1 0.00003920 0.0012395 -317.08
## + vs 1 0.00002130 0.0012574 -316.62
## + wt:hp 1 0.00000850 0.0012702 -316.30
## + drat 1 0.00000187 0.0012769 -316.13
## + am 1 0.00000185 0.0012769 -316.13
## + carb 1 0.00000000 0.0012787 -316.08
## + factor(gear) 2 0.00003001 0.0012487 -314.84
## + factor(cyl) 2 0.00001654 0.0012622 -314.50
## - hp 1 0.00046144 0.0017402 -310.22
## - wt 1 0.00221606 0.0034948 -287.91summary(car.step)##
## Call:
## lm(formula = I(1/mpg) ~ wt + hp, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0168690 -0.0049601 0.0007663 0.0040027 0.0142186
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.305e-03 4.094e-03 1.540 0.13435
## wt 1.149e-02 1.620e-03 7.089 8.45e-08 ***
## hp 7.479e-05 2.312e-05 3.235 0.00303 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.00664 on 29 degrees of freedom
## Multiple R-squared: 0.8471, Adjusted R-squared: 0.8365
## F-statistic: 80.33 on 2 and 29 DF, p-value: 1.494e-12Se obervă că se obține un model mult mai simplu, având ca termeni \(wt\) și \(hp\). Acesta este: \[1/pmg \sim wt + hp\]
Exercițiul 22.2
a.
#library("boot")
nuc.0 <- lm(cost ~ date + cap + pt + ne, data = nuclear)
summary(nuc.0)##
## Call:
## lm(formula = cost ~ date + cap + pt + ne, data = nuclear)
##
## Residuals:
## Min 1Q Median 3Q Max
## -157.894 -38.424 -2.493 35.363 267.445
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.756e+03 1.286e+03 -3.699 0.000975 ***
## date 7.102e+01 1.867e+01 3.804 0.000741 ***
## cap 4.198e-01 8.616e-02 4.873 4.29e-05 ***
## pt -1.289e+02 4.950e+01 -2.605 0.014761 *
## ne 9.940e+01 3.864e+01 2.573 0.015908 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 90.8 on 27 degrees of freedom
## Multiple R-squared: 0.7519, Adjusted R-squared: 0.7151
## F-statistic: 20.45 on 4 and 27 DF, p-value: 7.507e-08b.
plot(nuc.0, which = 1)
plot(nuc.0, which = 2)
Primul grafic trebuie analizat si verificat dacă reziduurile prezintă un trend de heteroscedasticitate. La o analizare atentă, pare că observațiile confirmă un trend de homoscedasticitate (sunt dispuse relativ aleator în jurul valorii 0). Observația \(19\) pare să nu se comformeze cu restul datelor, având o ușoară deviație, fapt confirmat și cel de-al doilea grafic unde se observă mai bine această deviere.
c.
size <- nrow(nuclear)
treshold <- 4 / size
plot(nuc.0, which = 4)
abline(h = treshold, lty = 3, col = 'red')
Se observă că există o observație care depășește pragul setat pentru distanța Cook (\(D_{i} > 4/n\)), aceasta fiind obervația \(19\). De la subpunctul precedent s-a constatat că aceasta nu se conformează cu dreapta de regresie liniară, iar prin metrica distanței Cook, acest lucru a fost confirmat ca fiind o observație cu un grad ridicat de influență asupra modelului de regresie.
d.
plot(nuc.0, which = 5, cook.levels = c(treshold, 0.5, 1))
Observația \(19\) este singura care se află în zona cu un grad de influență ridicat datorită valorii reziduale mari pe care o are. Observația \(10\) se află în vecinătatea acestei zone, delimitată de valoarea de prag \(4/n\), dar cum efectul de levier (leverage) nu este suficient de mare, deci aceasta nu are un grad mare de influență asupra modelului. Același lucru se poate spune și despre observația \(12\): dacă ar fi avut un efect de levier mai mare, atunci ar fi pătruns în zona punctelor cu influență ridicată.
e.
Observația \(19\) prezintă cea mai mare influență asupra modelului, deci ea va fi eliminată.
nuc.1 <- lm(cost ~ date + cap + pt + ne, data = nuclear[-19,])
summary(nuc.0)##
## Call:
## lm(formula = cost ~ date + cap + pt + ne, data = nuclear)
##
## Residuals:
## Min 1Q Median 3Q Max
## -157.894 -38.424 -2.493 35.363 267.445
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.756e+03 1.286e+03 -3.699 0.000975 ***
## date 7.102e+01 1.867e+01 3.804 0.000741 ***
## cap 4.198e-01 8.616e-02 4.873 4.29e-05 ***
## pt -1.289e+02 4.950e+01 -2.605 0.014761 *
## ne 9.940e+01 3.864e+01 2.573 0.015908 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 90.8 on 27 degrees of freedom
## Multiple R-squared: 0.7519, Adjusted R-squared: 0.7151
## F-statistic: 20.45 on 4 and 27 DF, p-value: 7.507e-08summary(nuc.1)##
## Call:
## lm(formula = cost ~ date + cap + pt + ne, data = nuclear[-19,
## ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -157.226 -46.025 -0.714 42.392 157.996
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.466e+03 1.046e+03 -4.270 0.000231 ***
## date 6.742e+01 1.518e+01 4.442 0.000146 ***
## cap 3.478e-01 7.235e-02 4.807 5.6e-05 ***
## pt -1.166e+02 4.029e+01 -2.894 0.007598 **
## ne 1.158e+02 3.164e+01 3.660 0.001128 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 73.68 on 26 degrees of freedom
## Multiple R-squared: 0.8027, Adjusted R-squared: 0.7723
## F-statistic: 26.44 on 4 and 26 DF, p-value: 7.868e-09plot(nuc.1,which=1)
plot(nuc.1,which=2)
După ce s-a eliminat observația \(19\), graficele reziduurilor s-au îmbunătățit în ceea ce privește îndeplinirea ipotezelor componentelor de eroare. Graficul **Residuals vs Fitted* prezintă o omogenitate mai mare în privința distribuției în jurul valorii 0. Pentru graficul \(Normal Q-Q\), valorile s-au apropiat de linia centrală. Se constată o îmbunătățire erorii reziduale standard: \(73.68\) vs \(90.8\) și pentru \(p-value\): \(7.868*10^{-9}\) vs \(7.507*10^{-8}\)
f.
#library("faraway")
#data("diabetes")
diab.0 <- lm(chol ~ age * frame + waist, data = diabetes)
summary(diab.0)##
## Call:
## lm(formula = chol ~ age * frame + waist, data = diabetes)
##
## Residuals:
## Min 1Q Median 3Q Max
## -125.156 -25.290 -4.328 21.261 221.876
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 128.7526 18.5098 6.956 1.53e-11 ***
## age 0.9555 0.2684 3.560 0.000418 ***
## framemedium 26.9255 15.5990 1.726 0.085139 .
## framelarge 35.6819 19.5654 1.824 0.068977 .
## waist 0.8343 0.4311 1.935 0.053706 .
## age:framemedium -0.3757 0.3364 -1.117 0.264839
## age:framelarge -0.7885 0.3784 -2.084 0.037844 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 42.22 on 381 degrees of freedom
## (15 observations deleted due to missingness)
## Multiple R-squared: 0.08764, Adjusted R-squared: 0.07328
## F-statistic: 6.1 on 6 and 381 DF, p-value: 4.068e-06Sunt \(15\) observații care au fost șterse: “(15 observations deleted due to missingness)”
g.
plot(diab.0, which = 1)
plot(diab.0, which = 2)
Primul grafic trebuie analizat si verificat dacă reziduurile prezintă un trend de heteroscedasticitate. La o analizare atentă, pare că observațiile confirmă un trend de homoscedasticitate (sunt dispuse relativ aleator în jurul valorii 0). Observațiile \(63\) și \(295\) par să aibă valori extreme, prezentând reziduuri mai mari. Acest lucru este confirmat și de al doilea grafic. Se observă ca valorile se află în apropierea dreptei, cu excepția unei tendințe de a se depărta în extremitatea dreapta. Observațiile \(63\) și \(295\) prezintă valori reziduale standardizate mari pentru care devierea de la dreaptă este evidentă.
h.
size <- nrow(diabetes) - 15
treshold <- 4/size
plot(diab.0, which = 5, cook.levels = c(treshold, 3 * treshold, 5 * treshold))
Valoarea de prag este \(4/388 = 0.0103\). Se observă că există mai multe puncte pentru care se află în zona observațiilor cu grad mare de influență. Există 3 puncte care depășesc de 3 ori valoarea prag aleasă, iar altele 3 în imediată vecinătate a acestei zone. Observațiile \(4\) și \(148\) au aproximativ același reziduu, însă având un efect de levier(leverage) ridicat, acestea pot fi categorisite ca fiind puternic influente. Observația \(63\) nu este un punct de pârghie puternic, dar având o valoare reziduală mare, se poate spune că aceasta are o influență ridicată asupra regresiei liniare.
j.
#exista probleme cu linkul si nu se returneaza mereu dataframe-ul. Trebuie comentate liniile 1808 -> ... in caz #ca se primeste eroare ca nu se poate citi raspunsul de la server
dia.url <- "http://www.amstat.org/publications/jse/v9n2/4Cdata.txt"
diamonds <- read.table(dia.url)
names(diamonds) <- c("Carat","Color","Clarity","Cert","Price")
diamonds.0 <- lm(Price ~ Carat + Color + Clarity + Cert, data = diamonds)
summary(diamonds.0)##
## Call:
## lm(formula = Price ~ Carat + Color + Clarity + Cert, data = diamonds)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1740.0 -428.8 -128.3 314.3 3634.1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 169.18 255.02 0.663 0.508
## Carat 12766.40 190.02 67.183 < 2e-16 ***
## ColorE -1439.09 207.98 -6.919 2.83e-11 ***
## ColorF -1841.69 195.23 -9.433 < 2e-16 ***
## ColorG -2176.67 200.39 -10.862 < 2e-16 ***
## ColorH -2747.15 202.91 -13.538 < 2e-16 ***
## ColorI -3313.10 212.71 -15.575 < 2e-16 ***
## ClarityVS1 -1474.57 159.67 -9.235 < 2e-16 ***
## ClarityVS2 -1792.01 171.19 -10.468 < 2e-16 ***
## ClarityVVS1 -689.29 159.93 -4.310 2.23e-05 ***
## ClarityVVS2 -1191.16 148.76 -8.007 2.73e-14 ***
## CertHRD 15.23 107.25 0.142 0.887
## CertIGI 141.26 128.26 1.101 0.272
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 710.4 on 295 degrees of freedom
## Multiple R-squared: 0.9581, Adjusted R-squared: 0.9564
## F-statistic: 562.5 on 12 and 295 DF, p-value: < 2.2e-16plot(diamonds.0, which = 1)
plot(diamonds.0, which = 2)
plot(diamonds.0, which = 3)
Pentru primul grafic Residuals vs Fitted se obervă că modelul nu este unul liniar. În al doilea grafic Normal Q-Q se observă devieri în partea dreaptă a modelului, culminând cu valori extreme pentru observațiile \(131\), \(116\) și \(279\). Normalitatea distribuției este afectată. În al treilea grafic se observă o oarecare deviație și repartizare constantă.
k.
diamonds.1 <- lm(log(Price) ~ Carat + Color + Clarity + Cert, data = diamonds)
summary(diamonds.1)##
## Call:
## lm(formula = log(Price) ~ Carat + Color + Clarity + Cert, data = diamonds)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.31236 -0.11520 0.01613 0.10833 0.36339
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.8011915 0.0496111 137.090 < 2e-16 ***
## Carat 2.8550130 0.0369676 77.230 < 2e-16 ***
## ColorE -0.0295093 0.0404611 -0.729 0.46638
## ColorF -0.1063582 0.0379807 -2.800 0.00544 **
## ColorG -0.2063494 0.0389848 -5.293 2.35e-07 ***
## ColorH -0.2878756 0.0394752 -7.293 2.81e-12 ***
## ColorI -0.4165565 0.0413818 -10.066 < 2e-16 ***
## ClarityVS1 -0.2019313 0.0310634 -6.501 3.41e-10 ***
## ClarityVS2 -0.2985406 0.0333027 -8.964 < 2e-16 ***
## ClarityVVS1 -0.0007058 0.0311121 -0.023 0.98192
## ClarityVVS2 -0.0966174 0.0289396 -3.339 0.00095 ***
## CertHRD -0.0088557 0.0208641 -0.424 0.67155
## CertIGI -0.1827107 0.0249516 -7.323 2.33e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1382 on 295 degrees of freedom
## Multiple R-squared: 0.9723, Adjusted R-squared: 0.9712
## F-statistic: 863.6 on 12 and 295 DF, p-value: < 2.2e-16plot(diamonds.1, which = 1)
plot(diamonds.1, which = 2)
plot(diamonds.1, which = 3)
În ciuda aplicării funcției logaritm asupra variabilei \(Price\), tot nu s-a liniarizat modelul. Se observă ca se obțin alte observații valori de extrem (\(214\), \(211\), \(152\))
l.
diamonds.2 <- lm(log(Price) ~ Carat + Color+ Clarity + Cert + I(Carat^2), data = diamonds)
summary(diamonds.2)##
## Call:
## lm(formula = log(Price) ~ Carat + Color + Clarity + Cert + I(Carat^2),
## data = diamonds)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.15411 -0.04120 -0.00911 0.04543 0.14158
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.075351 0.029202 208.049 < 2e-16 ***
## Carat 5.670616 0.079284 71.523 < 2e-16 ***
## ColorE -0.079247 0.017387 -4.558 7.59e-06 ***
## ColorF -0.155991 0.016328 -9.554 < 2e-16 ***
## ColorG -0.245033 0.016734 -14.643 < 2e-16 ***
## ColorH -0.339098 0.016969 -19.983 < 2e-16 ***
## ColorI -0.442606 0.017742 -24.947 < 2e-16 ***
## ClarityVS1 -0.244470 0.013358 -18.301 < 2e-16 ***
## ClarityVS2 -0.320183 0.014279 -22.424 < 2e-16 ***
## ClarityVVS1 -0.094008 0.013574 -6.926 2.74e-11 ***
## ClarityVVS2 -0.176701 0.012592 -14.032 < 2e-16 ***
## CertHRD -0.006223 0.008938 -0.696 0.4868
## CertIGI -0.025413 0.011536 -2.203 0.0284 *
## I(Carat^2) -2.102922 0.058022 -36.243 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0592 on 294 degrees of freedom
## Multiple R-squared: 0.9949, Adjusted R-squared: 0.9947
## F-statistic: 4445 on 13 and 294 DF, p-value: < 2.2e-16plot(diamonds.2,which=1)
plot(diamonds.2,which=2)
plot(diamonds.2,which=3)
Prin adăugarea termenului \(Carat^{2}\) în formulă se observă efectul de netezire a liniei de regresie. Acest al treilea model realizat este cea mai bună variantă dintre toate celelate.